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A  new  method  for  studying  electron-nuclear  dynamics  in  chemical  processes  is 
presented.  The  method  is  founded  on  the  Time-Dependent  Variational  Principle  and  the 
description  of  the  electronic  wavefuntions  via  coherent  states.  The  resulting  equations 
have  a  symplectic  structure  resembling  that  of  classical  mechanics. 

A  first  model  is  developed  with  the  nuclei  treated  in  the  limit  of  narrow  Gaussian 
wavepackets  and  the  electrons  restricted  to  a  single  determinantal  description.  The 
equations  of  motion  are  analyzed  to  understand  their  physical  significance.  These 
equations  are  used  to  study  the  p+H  and  a+He  collisions.  The  range  of  validity  of 
the  model  is  examined  and  future  developments  discussed. 
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CHAPTER  1 
INTRODUCTION 


Since  the  dawn  of  quantum  mechanics,  the  greatest  emphasis  has  been  placed  on  the 
solution  of  the  time-independent  Schrodinger  equation.  The  time-dependent  Schrddinger 
equation,  except  in  the  simplest  model  systems,  is  far  more  difficult  to  solve,  enhancing 
the  bias  toward  the  time-independent  studies.  Much  interesting  knowledge  has  been 
obtained  from  the  work  done  in  the  time-independent  picture,  and  more  will  come.  Yet 
a  time-dependent  scheme  gives  important  insight  on  the  actual  dynamics  and  specific 
mechanisms  important  for  the  evolution  of  a  physical  system.  Recent  advances  in 
computational  power  have  led  to  a  growing  interest  in  solving  time-dependent  problems. 
Among  the  phenomena  which  can  be  studied  with  a  time-dependent  method  are  molecular 
photodissociation,  collision  induced  dissociation  and  gas  phase  scattering.  Electron 
transfer  can  be  analyzed  in  detail  and  reaction  rates  can  be  computed.  Using  spectral 
analysis,  much  information  can  be  gleaned  from  a  time-dependent  wavefunction.  At  the 
same  time,  experiments  are  beginning  to  probe  chemical  reactions  at  the  femtosecond 
time  scale  [1],  which  no  time-independent  method  can  hope  to  describe.  Here  I  present  a 
theoretical  method  to  study  the  time  evolution  of  a  molecular  system  and  the  first  results 
obtained  from  its  use. 

The  potentials  considered  here  are  only  the  Coulombic  ones.  This  restriction  can 
easily  be  lifted  to  include  other  couplings  if  necessary.  Electrons  are  treated  quantum- 
mechanically  while  nuclei  can  be  treated  in  a  similar  fashion  or  semiclassically.  It  is  also 
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possible  to  treat  some  of  the  nuclei  classically  and  others  quantum  mechanically.  This 
can  be  important  for  such  processes  as  proton  transfer  through  quantum  tunneling  for 
which  it  is  not  necessary  to  have  a  fully  quantum  description  of  all  nuclei. 
The  nonrelativistic  molecular  Hamiltonian  is 


where  Hartree  atomic  units  have  been  used;  N  is  the  number  of  electrons  and  M  is  the 
number  of  nuclei.  The  nuclei  will  be  considered  as  point  particles.  The  equation  to  be 
solved  approximately  is  the  time-dependent  SchrOdinger  equation: 


In  most  cases,  attention  is  focused  on  doing  quantum,  semiclassical  or  classical 
nuclear  dynamics  on  electronic  potential  energy  surfaces.  This  involves  the  work  of 
up  to  three  research  groups  working  over  a  period  of  years.  One  group  determines  the 
potential  energy  surface  at  a  number  of  discrete  points.  Another  interpolates  the  results 
to  get  a  full  surface.  Finally  another  group  does  the  dynamics  on  the  surface.  However, 
it  is  possible  to  think  in  terms  of  a  more  general  approach,  to  develop  a  method  which 
in  principle  can  be  applied  to  any  molecular  system,  regardless  of  the  kind  of  properties 
that  are  to  be  analyzed,  and  in  which  the  three  steps  of  the  process  just  described  are 
done  in  one  step,  with  the  additional  bonus  that  the  electron-nuclear  dynamics  is  retained 


N   t  M      ,  N    M    ~         N    N  . 

,=1  1  A=l  LMA  i=l  A=l  r'A      .=1  »,  r" 


(l.D 


4w = w> 


(1.2) 


3 

fully.  With  such  a  method  it  is  possible  to  study  scattering  problems,  charge  transfer, 
follow  chemical  reactions,  predict  and  analyze  vibrational  spectra,  Raman  spectra,  etc. 
in  one  unified  manner. 

The  work  I  present  here  is  a  new  approach  to  such  a  general,  electron-nuclear 
dynamics  (END)  method.  The  method  has  been  developed  over  the  past  five  years.  It 
began  as  a  simple  model  of  an  electron  transfer  process  with  classical  nuclei  moving  on 
quartic  surfaces,  studied  by  Deumens,  Ohrn  and  Lathouwers  [2].  Then  I  did  some  work 
with  Deumens  and  Ohrn  to  describe  the  nuclei  quantum  mechanically  on  such  surfaces 
[3].  After  this  preliminary  work,  we  generalized  the  method  so  that  both  electrons  and 
nuclei  could  be  treated  dynamically.  This  led  to  a  computer  code  called  DYNAMO  from 
which  results  of  calculations  on  a  water  molecule  and  the  collision  process  of  a  proton 
on  a  hydrogen  atom  have  just  been  published  [4]. 

The  method  is  based  on  mathematical  methods  developed  over  the  past  30  years 
in  other  branches  of  physics.  Chapter  2  reviews  prior  time-dependent  schemes.  The 
theory  of  the  method  I  use  is  discussed  in  detail  in  Chapter  3  along  with  the  physical 
interpretation  of  the  resulting  equations.  Chapter  4  discusses  the  results  of  this  method 
when  applied  to  two  systems  I  have  studied.  The  first  is  a  proton  colliding  on  a  hydrogen 
atom,  for  collision  energies  varying  over  6  orders  of  magnitude  and  including  an  analysis 
of  the  electron  transfer  and  the  state-to-state  transfer  and  excitation  probabilities.  The 
second  is  the  scattering  of  an  alpha  particle  off  a  helium  atom  at  low  keV  energies. 


CHAPTER  2 

AN  OVERVIEW  OF  TIME-DEPENDENT  METHODS 
Potential  Energy  Surfaces 

A  widely  used  procedure  is  the  Born-Oppenheimer  approximation  (BOA).  It  is 
based  on  the  fact  that  nuclei  are  heavier  than  electrons,  and  move  more  slowly  than 
electrons  under  most  circumstances.  Thus,  to  a  good  approximation,  one  can  think  of  the 
electrons  as  moving  among  fixed  nuclei;  the  electrons  have  enough  time  to  accommodate 
themselves  to  the  changing  nuclear  degrees  of  freedom. 

With  this  concept,  a  particular  expansion  of  the  system  wavefunction  is  used.  The 

Hamiltonian  is  separated  as 
M 

2MA 

(2.1) 
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A  complete  set  of  eigenstates  of  the  electronic  Hamiltonian  is  found  and  the  system 
wavefunction  is  expanded  in  terms  of  these: 

*({*},  {Ra})  =  Erf*^}!  {^})^nucl({^}) 

(2.2) 

^^({r,},^})  =  £,elec({^})^lec({r;},{^}) 
No  approximation  is  made  up  to  this  point  Only  a  particular  expansion  of  the  system 
wavefunction  is  chosen  which  I  will  refer  to  as  a  Born-Oppenheimer  (BO)  expansion  of 
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states.  Allowing  the  Hamiltonian  to  operate  on  this  expansion  of  the  system  wavefunction, 
and  projecting  on  a  BO  electronic  state 


where 
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Eq.  (2.3)  defines  the  action  of  the  Hamiltonian  on  the  nuclear  wavefunctions.  The 
operators  T1  and  T2  are  dynamic  couplings  which  are  the  result  of  the  nuclear  momentum 
operator  acting  on  the  electronic  BO  states.  These  couplings  can  be  written  as 

M 


(2.5) 


and  if  the  set  of  electronic  states  is  complete, 

j\ —  1 

The  nuclear  momentum  operator  is  replaced  by 

-iVA  -»  -«VA  +  (2.7) 


in  Eq.  (2.3)  to  yield 
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This  Hamiltonian  for  the  nuclear  wavefunctions  is  no  longer  Hermitian,  since 
^-i'V/  +  P,^  is  not  Hermitian.  This  derives  from  the  expansion  of  the  full  wave- 
function  in  terms  of  the  electronic  eigenstates  [5]. 

Since  the  BO  electronic  states  are  instantaneously  in  equilibrium  with  the  nuclear 

— * 

configuration,  inertial  effects  between  electrons  and  nuclei  are  described  through  P^. 
The  P,4  couplings  represent  the  change  of  the  electronic  basis  functions  with  the  nuclear 
configuration.  They  contain  rotations,  distortions,  polarizations,  change  of  character  of 
the  electronic  basis  functions  as  well  as  the  change  of  the  electronic  basis  functions  due 
to  simple  displacement  of  the  nuclei  [6]. 

The  BO  approximation  (BOA)  neglects  the  P^  matrices,  effectively  decoupling  the 
equations  for  the  nuclear  wavefunctions.  Each  component  of  the  nuclear  wavefunction 
feels  an  effective  potential  given  by 


called  Potential  Energy  Surface  (PES).  Another  possibility  is  to  keep  the  diagonal  terms 
of  the  vector  potential  couplings  in  the  definition  of  the  PES.  This  is  usually  referred  to 
as  the  adiabatic  approximation. 


(2.9) 


The  definition  of  the  PES  and  the  use  of  the  Born-Oppenheimer  or  adiabatic  approx- 
imation is  the  starting  point  for  most  molecular  dynamics  today. 

Additional  approximations  are  usually  made,  such  as  treating  the  nuclei  purely  as 
classical  particles  on  these  surfaces,  or  to  treat  the  nuclei  as  quantum  particles  only 
insofar  as  to  quantize  the  vibrational  or  rotational  degrees  of  freedom  of  the  nuclei.  In 
the  past  decade,  efforts  have  been  made  to  describe  the  nuclei  in  terms  of  wavepackets 
moving  along  these  PES. 

The  breakdown  of  the  BOA  occurs  when  the  assumption  behind  the  approximation 
is  no  longer  valid,  i.e.  when  nuclear  motion  cannot  be  considered  to  be  slower  than 
electronic  motion,  or  when  PES  become  degenerate.  This  requires  that  dynamic  couplings 
be  retained  to  some  degree.  Processes  in  which  this  happens  are  avoided  crossings  of 
PES  (if  adiabatic  potentials  are  used,  or  crossings  for  diabatic  surfaces)  where  a  change 
in  the  electronic  character  of  the  states  may  generate  important  dynamic  couplings  [7]. 
Another  case  is  when  nuclei  are  vibrating  and  rotating;  as  nuclei  get  closer,  the  moment 
of  inertia  decreases,  increasing  the  angular  velocity.  Under  these  conditions  the  electrons 
may  not  have  time  to  keep  up  with  the  nuclear  motion,  leading  to  predissociation. 

To  generate  a  PES  for  dynamics  is  not  simple.  The  work  and  time  invested  is  very 
large.  Sometimes  it  takes  the  work  of  three  research  groups  working  over  a  period  of 
years  to  produce  the  electronic  energies,  the  interpolated  PES  and  do  the  dynamics.  Even 
then,  discrepancies  between  experiments  and  theory  can  become  painfully  obvious,  as 
witnessed  in  the  discussion  in  the  journals  over  the  existence  of  certain  resonances  in 
the  H+H2  scattering  problem  [8,  9]. 
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Time-Dependent  Methods  on  PES 


The  methods  used  today  to  do  dynamics  on  surfaces  can  be  divided  into  classical 
trajectory,  semiclassical  and  "exact"  quantum  mechanical.  Couplings  between  surfaces 
and  other  non  Born-Oppenheimer  effects  are  sometimes  treated  as  perturbations. 

Classical  Trajectory  Schemes 

The  simplest  method  is  the  classical  trajectory  scheme  [10-13]  whereby  nuclei  are 
propagated  as  classical  particles  on  a  single  PES.  The  greatest  amount  of  work  is  finding 
the  surface.  The  dynamics  is  simple  to  implement  and  the  results  are  simple  to  interpret. 
Extensive  sampling  of  phase  space  is  required  for  a  full  final  state  resolution,  i.e.,  to  have 
reaction  cross  sections  and  rates. 

Classical  trajectory  methods  are  capable  of  producing  good  results.  They  are  used 
to  describe  isolated  encounters  as  well  as  processes  in  a  condensed  phase  by  using  the 
generalized  Langevin  equation  [14,  15]. 

The  shortcomings  of  this  approach  are  the  neglect  of  quantum  effects,  such  as  zero- 
point  energy  and  tunneling,  which  may  be  important.  The  study  of  systems  with  a  single 
potential  energy  surface  of  importance  is  also  a  severe  restriction.  However,  the  most 
severe  limitation,  which  is  shared  by  all  methods  that  use  PES,  is  the  time  it  takes  to 
obtain  such  a  surface.  The  fitting  procedure  is  also  crucial,  since  errors  in  the  fitting  can 
lead  to  large  quantitative  errors  in  the  dynamics.  The  limitations  are  more  evident  as  the 
degrees  of  freedom  of  the  system  increase.  The  computational  effort  of  finding  points  on 
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the  surface  and  then  fitting  those  points  to  some  analytic  function  becomes  unmanageable 
for  any  system  having  more  than  three  or  four  atoms. 

Systems  investigated  with  the  classical  trajectories  include  the  H+H2  collision  [16] 
and  charge  transfer  in  the  H2+  +  H2  collisions  [17]. 
Semiclassical  Schemes 

The  semiclassical  approaches  try  to  correct  for  the  deficiencies  of  the  classical  tra- 
jectories without  going  to  the  full-blown  quantum  mechanical  approach.  There  are  many 
different  implementations,  from  time-independent,  such  as  Miller's  5-matrix  approach 
[18,  19]  in  which  the  Feynman  path  integral  is  evaluated  on  the  classical  trajectory  to  find 
an  approximation  to  the  5-matrix,  to  Heller's  Gaussian  wavepacket  path  integral  formula- 
tion of  semiclassical  dynamics  [20-22].  In  Heller's  work,  the  wavepacket  is  assumed  to 
remain  Gaussian  throughout  the  evolution,  with  the  form  exp[— a(x  —  q)2  +  ip(x  —  q)+c] 
and  equations  are  derived  for  the  parameters  q,p,a,  and  c.  The  wavepackets  are  propa- 
gated on  a  potential  energy  surface.  To  realize  the  propagation,  the  potential  is  quadrat- 
ically  expanded  about  the  instantaneous  center  of  the  wavepacket.  The  relationship  be- 
tween the  5-matrix  approach  and  the  time-dependent  wavepackets  has  been  studied  by 
Heller  [23]. 

To  correct  for  the  divergence  of  semiclassical  wavefunctions  near  caustics  a  gener- 
alization was  developed  [24,  25]  which  is  similar  to  Klauder's  global,  uniform  semiclas- 
sical approximation  for  wavefunctions  [26]  proposed  earlier.  In  this  work,  the  Gaussian 
wavepackets  are  generalized  so  that  all  the  time-dependent  parameters  are  complex  and 
the  manifold  of  states  labelled  by  p  and  q  is  searched  to  find  the  most  important  states 
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which  lead  to  the  same  final  state.  Other  approaches  along  similar  lines  include  the  work 
of  McDonald  [27]  and  Littlejohn  [28]. 

The  eikonal  approach  to  semiclassical  time-dependent  mechanics  has  been  proposed 
by  Micha  [29]  and  compared  to  Heller's  work  [30,  31],  yielding  very  similar  results. 
Applications  of  these  methods  include  work  on  the  photodissociation  of  methyl  iodide 
[30,  31].  A  comparison  between  the  semiclassical  approximation  and  an  "exact"  approach 
indicates  that  the  two  methods  agree  closely  for  this  system[31]. 

Another  semiclassical  computation  uses  Gaussian  wavepackets  in  the  interaction 
picture  for  the  H+H2  collision  [32].  It  has  been  compared  to  the  classical  trajectory 
work  [16].  The  two  methods  give  similar  results,  with  some  differences  in  the  details. 

Exact  Quantum  Mechanical  Schemes 

The  so-called  exact  methods  do  use  approximations.  Their  name  derives  from  the 
fact  that  the  errors  can,  in  principle,  be  made  as  small  as  desired. 

There  are  two  parts  to  consider:  the  representation  of  the  wavefunction  and  the 
actual  time  propagation  algorithm.  The  wavefunction  representation  determines  how  the 
operators,  such  as  the  Hamiltonian,  are  to  be  evaluated. 

Wavefunctions  are  represented  in  one  of  two  ways,  by  expanding  in  a  basis  set  or 
discretization  on  a  grid  of  points.  The  discretization  of  space  depends  on  the  PES  and 
the  choice  of  the  coordinates.  All  of  these  methods  work  only  with  internal  coordinates 
of  the  molecular  system,  avoiding  the  calculation  of  global  rotations  or  translations 
of  the  system.  The  effort  put  into  finding  good  internal  coordinates  is  nontrivial  [6]. 
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Transformations  to  such  coordinates  can  be  quite  a  problem.  It  is  much  more  convenient 
to  do  the  calculations  in  Cartesian  coordinates,  including  global  translation  and  rotation. 
For  this,  however,  it  is  necessary  that  the  conservation  laws  related  to  these  quantities 
remain  valid  within  the  approximations  made. 

The  Pseudo  Spectral  Fourier  Approximation  [33]  uses  the  wavefunction  represented 
on  a  grid  in  coordinate  space  and  then  uses  the  Discrete  Fourier  Transform  to  obtain 
the  momentum  space  representation.  Fast  Fourier  Transform  codes  allow  for  efficient 
switching  from  coordinate  to  momentum  space.  This  permits  fast  evaluation  of  operators 
such  as  the  potential  and  kinetic  operators  of  the  Hamiltonian. 

The  most  used  propagation  methods  can  be  divided  into  four  categories:  Second-order 
differencing  (SOD),  the  split-operator  method  (SO),  the  short-time  iterative  Lanczos  (SIL) 
method  and  the  Chebychev  expansion  (CE)  method. 

For  the  SOD  [34],  the  Schrodinger  equation  is  solved  by  approximating  the  time 
derivative  by  a  second  order  difference.  Short  time  steps  are  used,  and  the  wavefunction 
at  *n+i  is  evaluated  as 

0(tn+1)  =  V(<»-i)-2tA<tfV(<»)  •  (2.10) 

This  method  is  correct  through  second  order  in  the  time  step. 

The  other  three  methods  use  different  approximations  to  the  time  evolution  operator: 


U(e)  =  exp(-iHe) 


(2.11) 
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In  the  SO  scheme  [35]  the  exponential  operator  is  decomposed  symmetrically,  to  obtain 
an  approximate  evolution  operator  of  the  form 

U(e)  w  exp(-i|A')exp(-ieV)exp(-t'|/i:)  (2.12) 

which  is  also  correct  through  second  order. 
The  SIL  propagation  formula  [36]  is 

U(e)  «exp[-ieA(#,V>(0))]  (2-13) 

where  U  is  a  matrix  operator  in  the  Krylov  subspace,  which  is  generated  by  the 
Hamiltonian  and  the  initial  wavefunction,  and  A  is  the  tridiagonal  Lanczos  matrix 
representing  the  Hamiltonian  in  the  Krylov  space.  The  Krylov  subspace  is  generated 
through  the  action  of  the  Hamiltonian  operator  on  the  initial  wavefunction.  Thus  the 
N  -  I  Krylov  subspace  is  generated  by  the  N  vectors 

Uj  =  Hjxl>(0)   .  (2.14) 

The  Lanczos  recurrence  generates  a  set  of  orthogonal  polynomials  within  the  subspace 
which  represent  a  finite  polynomial  approximation  to  the  operator.  Then  the  projected 
subspace  representation  of  the  Hamiltonian  operator  has  a  tridiagonal  form.  The  operator 
is  diagonalized  and  the  propagation  is  performed  with  the  diagonal  eigenvalue  matrix. 
The  length  of  time  dictates  the  size  of  the  Krylov  space  needed  for  a  predefined  accuracy. 
Short  time  steps  are  used  in  order  not  to  lose  the  advantage  of  the  Lanczos  reduction. 
The  CE  method  [37]  approximates  the  exact  propagator  by  the  expansion 

N 

U(t)  =  ^^(t)Tn(-iHR)  (2.15) 

n=l 
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where  the  Hamiltonian  is  renormalized  so  that  its  spectrum  coincides  with  the  domain  of 
the  Chebychev  polynomials  T„.  The  normalization  is  carried  out  by  dividing  by  the  range 
of  energies  permitted  in  the  evolution,  AE  =  Emax  —  £min,  and  shifting  the  energy  level 
so  that  all  the  eigenvalues  fall  between  -1  and  1.  The  coefficients  a„  are  proportional 
to  Bessel  functions 

an(t)  =  2Jn(^lp)  (2.16) 

and  the  Chebychev  polynomials  Tn  are  obtained  from  the  recursion  relation 

to  =  V>(0) 

fa  =  -iHRil>(0)  (2.17) 

<f>n+l  =  -2iHR<f)n  +  <f>n-\  . 

The  CE  method  is  a  global  propagator.  It  is  such  that  the  number  of  terms  in  the 
expansion  does  not  decrease  significantly  for  small  time  steps.  Therefore,  for  efficiency, 
the  time  step  should  be  large,  sometimes  to  the  point  that  a  single  time  step  completes 
the  calculation.  Restrictions  of  this  method  are  that  only  time-independent  Hamiltonians 
may  be  used.  Differences  with  the  Lanczos  scheme  include  the  fact  that  the  coefficients 
are  fixed  for  the  CE  while  in  Lanczos  they  depend  on  the  initial  state.  Also,  the  SIL  uses 
small  time  steps  for  efficiency  compared  to  the  long  time  steps  of  the  CE  method. 

Specific  advantages  of  each  of  these  methods  are  reviewed  by  Kosloff  [38].  One 
of  the  limitations  for  all  these  methods  is  the  number  of  atoms  in  a  molecular  system 
that  can  be  studied  due  to  computational  complexities.  At  this  point  systems  of  up  to  3 
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atoms  are  possible,  with  some  work  being  started  on  systems  of  4  atoms  in  which  some 
internal  degrees  of  freedom  can  be  frozen. 

A  few  representative  works  with  "exact"  methods  are  reactions  [39]  and  absorption 
spectra  [40].  Because  of  the  SOD's  ease  of  implementation  and  satisfactory  accuracy  it 
has  been  used  extensively  for  a  variety  of  problems,  including  eigenspectra  [41],  non- 
adiabatically  coupled  systems  [42],  photodetachment  spectra  [43]  and  systems  with  time- 
dependent  Hamiltonians  [44].  The  CE  method  has  been  used  on  atom -diatom  collisions 
[45],  photodissociation  [46]  and  computation  of  energy  levels  [47,  48]. 

To  reduce  the  number  of  grid  points  that  are  required  for  the  calculations,  most  work 
is  done  with  absorbing  potentials  which  are  placed  at  the  boundary  of  the  coordinate 
space  considered  to  be  important  for  the  dynamics  of  the  reaction. 

Time-Dependent  Methods  Without  PES 


Most  chemistry  does  not  happen  on  a  single  potential  energy  surface  with  only  three 
or  four  atoms.  Determining  PES  for  higher  numbers  of  atoms  becomes  a  difficult,  if  not 
impossible  task.  These  facts  make  the  search  for  a  method  that  does  not  use  surfaces 
more  urgent. 

The  limitations  of  PES  have  generated  interest  in  developing  methods  that  avoid 
their  use.  One  example  in  which  methods  using  PES  fail  is  the  collision  of  a  proton  with 
a  hydrogen  atom  at  energies  above  10  eV.  In  this  case  many  PES  would  be  required, 
since  excitation  to  n=2  states  becomes  important  and  provides  for  interesting  physical 
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phenomena.  The  collision  has  been  studied  on  the  lowest  energy  surfaces  made  up  by 
the  molecular  ls<r  and  2pa  states  [49-51],  but  for  energies  below  10  eV. 
Close-Coupling  Methods 

For  collisions  at  high  (above  1  keV)  energies,  the  close-coupling  methods  developed 
by  Fritsch,  Lin,  Kimura,  Thorson  and  Lane  and  others  [52-56]  have  been  successful 
in  describing  collisions  for  systems  with  one  electron,  or  in  which  all  but  one  electron 
can  be  approximated  by  some  potential.  Work  on  systems  which  explicitly  retain  two 
electrons  is  new  [57]  and  very  time  consuming.  The  close-coupling  method  projects  the 
time-dependent  Schrodinger  equation  on  a  basis.  The  nuclei  are  not  treated  dynamically 
but  trajectories  are  prescribed  for  them.  These  are  usually  either  straight  line  or  Coulomb 
trajectories.  The  nuclear  motion  with  this  approximation  leads  to  a  time-dependent 
potential  for  the  electrons.  It  is  only  the  expansion  coefficients  of  the  electronic  states 
that  are  calculated  for  different  times. 

Atomic  or  molecular  orbitals  with  electron  translation  factors  (see  the  end  of  the 
chapter  for  a  description  of  these  factors)  are  used  as  the  basis.  For  more  than  one  electron, 
the  basis  is  made  up  of  different  determinants  of  atomic  or  molecular  orbitals.  The 
electronic  states  are  described  as  a  CI  wavefunction  with  time -dependent  CI  coefficients. 
For  this  reason  as  the  number  of  electrons  grows,  the  method  becomes  more  and  more 
difficult  to  use. 

There  are  several  implementations,  such  as  Fritsch  and  Lin's  AO+  [54,  55]  or  the 
AO-MO  matching  procedure  of  Kimura  and  Lin  and  Winter  and  Lane  [58,  59].  The 
differences  depend  on  the  basis  sets  being  used.  In  general  the  electronic  wavefunction 
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is  written  as 


N 


!*(<)>  =  Efl*(')W)) 


(2.18) 


for  a  basis  set  of  N  functions.  This  basis  set  changes  with  time,  since  the  basis  functions 
follow  the  nuclei.  Projecting  the  Schrodinger  equation  on  the  truncated  basis  gives 


where  the  prescribed  nuclear  trajectories  make  the  electron-nuclear  attraction  and  nuclear- 
nuclear  repulsion  potentials  time-dependent.  The  equations  take  the  form  [53] 


where  the  final  expression  for  N  and  M  will  depend  on  the  basis  chosen  and  the  prescribed 
trajectory. 

Since  the  close-coupling  method  does  not  account  for  nuclear  dynamics  it  cannot  be 
successfully  used  for  collision  energies  below  about  1  keV.  The  path  followed  by  the 
nuclei  is  very  important  in  determining  transitions  at  lower  energies[60].  For  energies 
below  1  keV  it  is  essential  to  include  electron-nuclear  dynamics.  The  close-coupling 
methods  have  never  been  used  for  systems  in  which  more  than  two  electrons  are  explicitly 
described. 

Methods  with  Nuclear  Dynamics 

Several  methods  have  been  proposed  in  this  field.  One  method,  the  Car-Parrinello 
approach  [61],  has  gained  much  attention  among  theoretical  chemists  because  of  its 


(v*(oi4  -  #i*> = ° 


(2.19) 


(2.20) 


Njk{i)  =  tohfo),   Mjk  =  Wjlij-  -  H(t)\^k) 
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applicability  to  larger  systems.  The  method  uses  an  optimization  technique  to  follow  the 
ground  state  surface  and  uses  the  Hellman-Feynman  theorem  to  compute  the  forces  on 
the  nuclei.  It  is  an  approximation  to  dynamics  on  a  surface  which  requires  calculating 
only  the  part  of  the  PES  needed.  There  is  no  firm  theoretical  work  which  shows  fully 
why  the  method  works,  in  particular  why  it  works  better  than  if  an  SCF  calculation  were 
to  be  performed  at  each  position  instead  of  using  a  molecular  dynamics  optimization 
procedure  [62].  Hartke  and  Carter  use  this  method  on  Na4  singlet  and  triplet  states  [63]. 

Field  [64]  has  used  the  Time-Dependent  Hartree-Fock  method  for  a  closed-shell 
restricted  Hartree-Fock  wavefunctions  using  the  neglect  of  diatomic  differential  overlap 
(NDDO)  approximation.  In  his  approach,  the  TDHF  equations  drive  the  electronic 
degrees  of  freedom,  subject  to  the  constraint  of  normalized  molecular  orbitals: 


where  the  integration  over  spin  has  been  performed,  c  are  the  expansion  coefficients  of 
the  molecular  orbitals  in  terms  of  an  orthogonal  basis  and  F  is  the  Fock  matrix  in  the 
basis.  The  nuclei  are  assumed  to  evolve  as  in  classical  dynamics,  with  the  forces  given 
by  the  gradient  of  the  energy  of  the  system 


(2.21) 


Fji  =  c]Fc, 


dRj_ 
dt 


— * 

Mi 


(2.22) 
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The  energy  E  is  expressed  as 


(2.23) 

#ekc  =  22^cJi<ViJV+   ^2   c*fiic>>ic*<rjcXj(2(fl(T\1/X)  ~  (/H-W) 

Field  has  applied  this  method  to  the  study  of  properties  such  as  dipole  moments  and 
power  spectra  for  lithium  hydride,  water  and  formaldehyde  [64]. 

In  Field's  approach  discussed  up  to  this  point,  there  is  a  curious  omission.  The 
electronic  equations  of  motion  are  derived  without  considering  the  translation  of  the 
basis  set  with  the  nuclei.  This  neglects  couplings  in  the  equations,  couplings  which 
transform  the  electronic  coefficients  between  basis  sets  at  different  times. 

Another  approach  has  been  developed  by  Runge,  Micha  and  Feng[60].  Starting  from 
the  TDHF  approximation,  and  making  no  further  approximations  such  as  the  NDDO  or 
restricted  close-shell  determinant  of  Field,  they  derive  equations  of  motion  for  a  basis 
which  includes  electron  translation  factors.  In  this  case  the  transformation  due  to  the 
translation  of  the  basis  is  retained  initially.  However,  by  using  electron  translation  factors 
it  is  eliminated  because  the  coefficients  are  written  in  terms  of  a  moving  basis.  They 
do  neglect  another  coupling,  namely  the  one  between  the  nuclear  accelerations  and  the 
electronic  degrees  of  freedom,  which  derives  from  the  electron  translation  factors.  In  their 
first  implementation  (Ref.  60)  they  also  neglect  the  electron  translation  factors  in  the 
actual  calculation  of  the  Fock  matrix.  Thus,  electron  translation  factors  are  used  only  to 
eliminate  the  basis  transformation  effects.  Because  of  these  approximations  the  equations 


19 

reduce  to  ones  similar  to  those  of  Field  but  without  the  additional  approximations  of 
NDDO  and  restricted  close- shell  determinant.  Runge  et  al  solve  an  equation  for  the 
density  matrix  instead  of  the  orbital  coefficients.  The  equations  to  be  solved  are  [60] 

iP^  =  S-1F7P7-P7F7S1 

i 

(2-24) 

dRj_  Pj_ 
dt  ~  Mi 

dP]  dE 
dt  ~  &R7 

where  E  is  the  total  energy  of  the  system.  This  method  has  been  applied  to  systems  with 
one  electron,  such  as  the  p+H  collision  [60].  Work  in  progress  by  this  group  lifts  the 
neglect  of  electron  translation  factors  in  the  calculation  of  the  Fock  matrix. 

Another  approach  is  taken  by  Meyer  and  Miller  [65].  They  too  derive  equations 
which  couple  the  nuclear  and  electronic  degrees  of  freedom.  Their  approach  is  formal 
and  has  not  been  used  in  an  ab  initio  fashion  on  any  real  systems.  They  assume  the 
existence  of  a  complete  diabatic  basis  for  the  electrons  (i.e.  one  which  does  not  depend 
on  the  nuclear  positions).  This  has  the  effect  of  decoupling  the  electronic  and  nuclear 
degrees  of  freedom,  except  through  the  electron-nuclear  attraction.  The  basis  considered 
is  not  made  up  from  an  expansion  of  one-electron  orbitals,  but  is  a  general  TV-electron 
basis.  The  electronic  degrees  of  freedom  are  the  expansion  coefficients  of  the  electronic 
state  in  the  basis.  The  coefficients  are  written  as  a  norm  times  a  phase  factor.  The  nuclei 
are  approximated  by  point-like  classical  particles.  Using  the  time-dependent  variational 
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principle  they  derive  the  equations  of  motion.  The  equations  for  the  electronic  degrees 
of  freedom  look  just  like  the  action-angle  equations  of  motion  of  classical  mechanics. 
Thus  the  whole  system  of  equations  looks  very  "classical",  although  the  electrons  are 
treated  quantum-mechanically.  A  canonical  transformation  can  be  performed  to  obtain 
the  equations  of  motion  for  the  electrons  and  nuclei  using  an  adiabatic  electronic  basis. 
In  this  case  the  equations  of  motion  for  the  electrons  and  nuclei  are  [65] 

dpi 

■  _V/v  (p  +  p) 

Pl~    dxi~     V       dx>      dx>  M« 


(2.25) 


•        dH  dP    (p  +  P) 

Qk  =  -KTT  =  EK  + 


NK  = 


dNK  dNK  Mi 

dH         dP    (p  +  P) 


dQK       dQK  M, 
with  the  wavefunction  written  as 

|0)  =  X^xpHQ*)^!*)  (2-26) 
K 

where  x,  p  are  the  classical  nuclear  position  and  momenta,  P  are  the  momentum  couplings 
defined  in  the  BO  Hamiltonian  (Eq.  2.8),  |A')  are  the  adiabatic  electronic  states  and  the 
Hamiltonian  has  been  taken  to  the  classical  limit  for  the  nuclear  degrees  of  freedom. 
Electron  Translation  Factors 

The  Perturbed  Stationary  State  method  (PSS)  [6,  66]  is  used  for  many  atomic  and 
molecular  collisions  calculations.  Starting  from  the  BO  expansion  of  the  Hamiltonian, 
Eq.  (2.8),  this  method  does  not  ignore  the  dynamic  couplings  like  the  BOA.  Delos  [6] 
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lists  shortcomings  of  the  PSS  approach  when  the  electronic  as  well  as  nuclear  basis  is 
truncated,  i.e.  when  approximations  are  made  to  the  BO  expansion.  He  notes  that  the 
expression  for  T2  in  terms  of  P  in  Eq.  (2.6)  is  not  valid  for  a  truncated  basis,  though 
he  states  that  it  is  usually  sufficiently  accurate  within  the  approximation.  More  serious 
problems  he  notes  are: 

1.  The  individual  terms  in  the  BO  expansion  of  the  system  state  do  not  satisfy  the 
scattering  boundary  conditions. 

2.  The  P  matrix  elements  contain  couplings  of  infinite  range. 

3.  P  matrices  contains  fictitious  "origin  dependent  "  couplings. 

4.  Matrix  elements  in  this  formulation  do  not  contain  momentum  transfer  factors  which 
are  needed  to  describe  Doppler  shifts  in  the  energy  spectrum  of  moving  systems. 

Delos  suggests  a  solution  via  the  so-called  electron  translation  factors  (ETF's),  which 
eliminate  many  of  these  problems.  The  solution  he  proposes  is  essentially  to  include  in 
each  electronic  basis  function  a  complex  phase  factor  which  cancels  some  of  the  offensive 
terms  listed  above.  This  is  effectively  a  correction  to  the  BO  expansion,  accounting  for 
the  momentum  of  the  electrons  in  the  process,  not  through  the  dynamic  couplings  but 
through  the  basis. 

This  correction  is  used  in  the  close-coupling  method  for  atomic  collisions,  and  also 
in  the  approach  by  Runge,  Micha  and  Feng  to  electron-nuclear  dynamics. 

The  following  example  should  clarify  the  issue.  If  an  H  atom  is  moving  with  respect 
to  an  inertial  frame  (which  I  will  call  the  lab  frame),  its  ground  state  should  satisfy  the 
time-dependent  SchrOdinger  equation.  It  is  a  simple  exercise  to  show  that  the  orbital 
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exp(ime(u  •  f-  v2t/2)  -  eit)<f>is(?—  Rj  where  v  is  the  velocity  of  the  moving  atom, 
R  its  instantaneous  position  and  <f>is  the  time-independent  ground  state  solution  with 
energy  t\  solves  the  problem.  To  get  a  reasonable  approximation  of  this  orbital  as  an 
expansion  in  terms  of  a  truncated  set  of  the  lab  frame  eigenstates  of  the  H  atom,  then 
orbitals  with  n  greater  than  1  are  necessary.  The  greater  the  velocity,  the  greater  the 
size  of  the  basis  needed.  If  the  basis  used  for  calculations  is  fixed  to  the  lab  (even  if  it 
"follows"  the  nuclear  positions,  i.e.  using  an  instantaneous  basis  centered  at  the  position 
of  the  nuclei,  but  which  is  not  moving  relative  to  the  lab  frame),  a  small  basis  will  not 
appropriately  describe  the  ground  state  of  the  moving  H  atom.  Another  way  of  viewing 
this  is  that  the  description  is  not  Galilean  invariant  when  using  a  truncated  basis  since 
calculations  done  with  the  basis  fixed  in  a  different  frame  gives  different  results. 

The  solution  to  the  problem  is  simple  for  atoms  [6].  One  just  uses  the  phase  factors 
as  shown  above  on  each  basis  function  associated  with  an  atom  and  Galilean  invariance 
is  recovered  even  for  a  truncated  basis.  If  the  calculations  use  molecular  orbitals  instead 
of  atomic  orbitals,  the  problem  becomes  a  difficult  one.  In  fact,  no  fully  satisfactory 
solution  has  yet  been  found  to  this  problem.  The  ETF's  in  this  case  have  to  be  such  that 
if  the  electron  is  close  to  a  given  nucleus  it  has  the  appropriate  velocity  for  that  nucleus, 
with  a  continuous  transformation  of  the  velocity  between  nuclear  centers.  Several  forms 
have  been  proposed  [56,  67],  but  each  one  has  its  drawbacks.  Another  problem  with  the 
MO-ETF's  is  that  calculations  depend  on  the  origin  of  the  electronic  coordinates.  After 
more  then  a  decade  of  working  on  this  problem,  groups  using  close-coupling  methods 
have  settled  on  the  AO-MO  matching  scheme  [57].  In  this  scheme  the  orbitals  used 
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outside  a  certain  distance  between  nuclei  are  atomic  and  contain  the  electron  translation 
factors.  At  the  limiting  distance  these  AO's  are  matched  to  molecular  orbitals  without 
electron  translation  factors  to  continue  the  calculations  with  the  PS  S.  When  the  fragments 
separate,  at  the  limiting  distance  the  matching  from  the  MO  to  the  AO  basis  is  again 
performed.  In  this  way  Galilean  invariance  is  preserved  outside  the  interaction  region. 

An  approach  where  the  MO's  are  described  in  terms  of  an  LCAO  expansion  with 
AO's  having  ETF's  was  attempted,  but  failed  (A.  Riera,  personal  communication).  The 
failure  stems  from  the  fact  that  the  Hamiltonian  carries  no  information  about  the  moving 
frame  and  so  the  MO's  obtained  undo  the  effects  of  the  ETF's. 

When  ETF's  are  used,  there  appears  to  be  no  work  which  has  included  the  cou- 
plings that  arise  for  accelerated  nuclei.  For  close-coupling  treatments  with  straight  line 
trajectories  and  constant  velocities  these  couplings  are  zero. 


CHAPTER  3 
THE  END  FORMALISM 


The  fully  quantum,  exact  methods  are  setting  benchmarks,  but  the  use  of  the  potential 
energy  surfaces  severely  restricts  them,  a  problem  shared  by  all  methods  relying  on  PES. 
Even  four  atom  systems  are  beyond  reach,  greatly  limiting  the  chemical  systems  that 
can  be  studied.  The  move  towards  the  elimination  of  potential  energy  surfaces  is  very 
recent.  As  already  discussed,  several  proposals  have  been  made.  However,  either  the 
theory  behind  the  methods  is  not  well  understood  or  the  approximations  made  are  many, 
and  at  times  confusing.  Here  I  present  another  method  which  does  not  rely  on  PES.  The 
approximations  made  are  few  and  simple  to  understand. 

The  chapter  is  divided  in  four  parts.  The  first  is  the  description  of  the  Time-Dependent 
Variational  Principle.  The  second  is  the  coherent  state  formalism.  In  the  third  part  these 
tools  are  used  to  find  the  equations  for  electron-nuclear  dynamics  (the  END  equations). 
There  is  also  a  description  of  the  approximations  made  in  the  first  model  developed.  The 
fourth  part  deals  with  the  tools  developed  to  analyze  the  results  of  a  calculation. 

The  Time-Dependent  Variational  Principle 

Three  main  time-dependent  variational  principles  are  used  today.  They  are  the 
McLachlan  variational  principle  [68],  the  Dirac-Frenkel  variational  principle  [69,  70]  and 
the  Time- Dependent  Variational  Principle  (TDVP)  [71].  All  of  these  have  been  shown  to 
be  equivalent  [72]  if  the  trial  wavefunction  is  described  with  complex  parameters  and  is 
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analytic  in  these  parameters.  Because  the  form  of  the  wavefunction  used  here  satisfies 
these  two  conditions,  the  three  variational  principles  will  give  the  same  results.  I  will 
describe  and  use  the  TDVP. 

If  we  assume  the  existence  of  a  set  of  electronic  and  nuclear  parameters  £  which 
completely  define  the  quantum  state  of  the  system,  I  will  label  the  state  with  these 
parameters:  \().  Using  Hartree  atomic  units  so  that  h  =  1,  the  following  Lagrangian 
is  proposed: 

=  (CIA'IO(CIC)-1 

where  the  left-acting  time  derivative  is  noted  as  ^,  and  the  operator  K  is  defined.  The 
action  is  then 

A  =  I*  Ldt 
Jti 

(3.2) 

and  the  TDVP  requires  that 

sa = 8C  [K(ci^ic>  _  ^((c|)ic>) _  (c|//|c>] (cic>_1^ = °  ■  (33) 

The  TDVP  will  result  in  the  Schrodinger  equation  if  the  trial  wavefunction  used  has  no 

restriction  in  the  Hilbert  space  of  the  system.  With  a  partial  integration: 

r<2 


SA  =  J^[(6C\t^\0-(SC\H\0 


(3.4) 


-W\Q{C\Q-l{6C\Q  +  c.c.  ((IQ-'dt 


i-i. 
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Because  the  variation  is  arbitrary 


10 


(3.5) 


which  is  the  Schrodinger  equation  if  the  right  hand  side  is  zero.  If  the  wavefunction  is 
multiplied  by  global  time-dependent  phase  factor  exp(z'7)  with  7  satisfying  the  equation 

*  "  IcicT  (  } 

then  it  is  simple  to  see  that  the  variation  of  the  action  results  in 


recovering  the  Schrodinger  equation.  The  global  phase  is  important  in  time -correlation 
functions. 

In  practical  applications  the  trial  wavefunction  is  restricted  to  a  submanifold  of  the 
full  Hilbert  space  and  an  approximation  to  the  time-dependent  equation  is  solved. 
Introducing  the  notations 


(3.7) 


=  exp  (t'7)|C) 


S(C',0  =  «I0 


(3.8) 


£(C,C)  =  <C|tf|0/(CI0 


then 


(3.9) 


This  differential  equation  can  be  integrated  with  a  simple  quadrature. 
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Because  all  time  dependence  of  the  state  is  through  the  parameters,  Eq.  (3.3)  is 
equivalent  to 


dt 


-f 


dt 


-/?[(-?S6-S)* 


(3.10) 


dt 


The  second  step  in  Eq.  (3.10)  involves  a  partial  integration  of  some  terms  with  respect 
to  t.  Since  all  the  8(  and  their  complex  conjugates  are  independent  variations  we  obtain 
the  dynamical  equations 


(3.11) 


with  the  elements  of  the  metric  matrix  defined  as 


dE 


d2\nS 


(3.12) 


The  matrix  C  is  clearly  Hermitian.  The  phase  factor  does  not  influence  the  evolution  of 
the  other  parameters,  and  can  be  computed  afterwards.  In  practice,  it  is  convenient  to 
integrate  (3.9)  along  with  the  equations  (3.11)  for  the  other  parameters. 
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The  equations  (3.11)  can  be  written  in  matrix  form  as 

<s  mnt> 

Furthermore,  we  can  define  a  Poisson  bracket  {,}  for  two  functions  /  and  g  depending 
on  (  and  (*  as 

(/.f}-<#  &>(T  .c'->)(f)-  (3-M) 

With  the  symplectic  structure  defined  by  this  Poisson  bracket,  the  evolution  equations 
(3.11)  assume  the  standard  form 

C  =  {C,£}, 

(3.15) 

which  shows  that  the  energy  E  assumes  the  role  of  Hamiltonian  or  generator  of  time 

translations,  while  (  and  C*  are  conjugate  variables, 

{C,C}  =  0, 

{C*,C*}=0,  (3.16) 

{c,n  =  -ic-1. 

The  phase  space  is  fiat  only  when  C  is  the  unit  matrix. 

Coherent  States 

The  parametrization  of  the  wavefunction  is  of  crucial  importance.  A  poor  parametriza- 
tion  will  lead  to  a  complicated  phase  space  with  a  metric  which  may  lead  to  integration 
problems  and  additionally  make  the  equations  difficult  to  understand.  This  problem  was 
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first  addressed  by  Thouless  [73]  for  nuclear  systems  using  the  TDHF.  The  main  problem 
he  faced  is  that  a  single  determinant  is  unchanged  under  a  unitary  transformation  of  its 
orbitals.  This  means  that  there  are  redundant  parameters  in  the  definition  of  a  deter- 
minantal  state.  The  parametrization  he  found  turns  out  to  be  a  special  case  of  a  more 
general  scheme,  referred  to  as  coherent  states. 

Schrodinger  first  proposed  the  concept  of  coherent  states  in  1926  [74]  in  connection 
with  the  classical  harmonic  oscillator.  Coherent  states  were  later  used  by  Glauber 
[75-77]  and  Sudarshan  [78]  in  quantum  optics  in  the  early  60's.  Since  then  there 
has  been  interest  in  them  from  both  a  mathematical  [79]  and  physical  [80]  point  of 
view.  The  representations  of  Lie  groups  can  be  studied  with  coherent  states,  and  the 
range  of  physical  systems  treated  with  this  technique  goes  from  condensed  matter  and 
thermodynamics  to  elementary  particle  physics  and  path-integral  developments,  as  well 
as  to  study  the  relationship  between  quantum  and  classical  dynamics  (for  a  compilation 
of  these  developments  see  Klauder  and  Skagerstam's  book,  Ref.  80).  Coherent  states  are 
also  useful  in  the  study  of  unrestricted  Hartree-Fock  theory  applied  to  molecules  [81]. 

A  coherent  state  satisfies  two  properties  [80]: 

1.  A  coherent  state  |^)  is  a  strongly  continuous  function  of  the  label  z.  This  means  that 

lim  I  \z)  -  \z')  I  =  0  (3.17) 

z— >z' 

2.  There  is  a  positive  measure  6z  on  the  space  Q,  of  labels  z  such  that 

/  =  /  \z)(z\6z  (3.18) 
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The  first  property  indicates  that  there  is  a  one-to-one  mapping  between  coherent  states 
and  points  on  the  label  space  fi.  The  second  property  shows  that  the  coherent  states  are 
a  complete  set,  since  any  state  can  be  decomposed  in  terms  of  coherent  states.  Because 
of  the  continuous  label,  the  set  cannot  be  an  orthogonal  basis,  but  must  be  overcomplete. 

Although  there  are  many  forms  of  coherent  states  [80],  I  will  restrict  myself  to 
ones  related  to  compact  Lie  groups.  Let  U(g)  be  a  continuous,  irreducible  unitary 
representation  of  a  compact  Lie  group  G  and  |0)  be  some  normalized  reference  state, 
also  called  the  extremal  state.  Then  the  states  generated  by  the  operation  of  the  group 
on  the  reference, 


satisfy  the  first  property  from  (g\g')  =  (0|C/(sr_1p')|0)  and  the  continuity  properties  of 
Lie  groups.  To  show  that  the  second  property  is  also  satisfied,  take 


where  dg  is  the  group  invariant  measure.  By  virtue  of  this  measure,  it  is  clear  that  for 
any  g',  PU(g')  =  U(g')P.  Schur's  Lemma  then  requires  that  P  be  proportional  to  the 
identity.  Renormalizing  the  invariant  measure  by  a  constant  leads  to  P  =  7. 

The  states  \g)  may  not  all  be  distinct.  To  obtain  a  distinct  states  let  h  be  an  element 
of  a  subgroup  H  of  G  defined  by  the  property  that  these  group  elements  only  change 
the  reference  by  a  phase  factor,  i.e. 


\g)  =  U{g)\0) 


(3.19) 


(3.20) 


U{h)\0)  =  exp[ia{h)}\0) 


(3.21) 
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H  is  called  the  stability  subgroup  and  the  coset  G/H  provides  a  unique  decomposition 
of  any  element  g  belonging  to  G  [82]: 

g  =  ch  (3.22) 

where  c  is  a  coset  representative  and  h  belongs  to  the  stability  subgroup.  Then 

U(g)\0)  =  U(c)U(h)\0) 

=  U(c)\0)exp[ia(h)]  (3.23) 
=  |c)  exp  [ic*(/i)] 

and  |c)  is  the  coherent  state,  defined  through  the  action  of  the  unique  coset  representative 

c. 

For  a  single  determinant,  the  stability  subgroup  is  the  subgroup  made  up  of  unitary 
transformations  which  only  mix  the  occupied  or  unoccupied  orbitals.  These  do  not  change 
the  determinant.  By  using  the  coset  parameter  representation  the  possibility  of  parameters 
changing  without  changing  the  determinant  is  eliminated  from  the  outset. 

I  will  return  to  the  single  determinant  a  little  further  on.  First  I  will  illustrate  how  to 
apply  the  above  strategy  to  the  simplest  coherent  state,  the  harmonic  oscillator  coherent 
state,  starting  from  the  Heisenberg-Weyl  Lie  algebra.  These  coherent  states,  also  called 
canonical  coherent  states  [80]  can  be  used  as  a  first  approximation  to  treating  nuclei 
quantum  mechanically  for  molecular  dynamics. 


32 

The  Heisenberg-Weyl  algebra  is  spanned  by  the  operators  /,  a,  al  and  n  =  a^a  with 
the  commutation  relations: 


n,  a1 


[n,  I]  =  0, 


[n,a]  =  -a, 


=  0, 


(3.24) 


a,  a 


=  /,        M  =  o. 

The  carrier  space  of  irreducible  representations  for  this  Heisenberg-Weyl  group  (H4)  is 

spanned  by  the  number  eigenstates  of  n: 

n\k)  =  k\k),     k  =  0,1,2,.... 

k  (3.25) 


jot)' 
(*!)' 


I*)  =  77^1°) 


The  state  |0)  will  be  used  as  the  reference.  It  is  the  lowest  weight  state  for  the  number 
operator  n.  The  group  elements  h  =  exp  [i(Sn  +  <f>I)]  of  #4  only  change  the  reference 
by  a  phase  factor  exp(i<f>).  These  elements  form  the  stability  subgroup  U{1)  <g>  U(l)  of 
#4.  A  general  element  of  H4  can  be  written  as 


g  =  exp  ^aa*  —  a*a^  exp  [i(8n  +  <j>I)] 


(3.26) 


so  that  the  coherent  state  is 

\a)  =  exp  (aa^  —  a*aj  |0) 


=  exp  ^— |a|2/2^  exp  ^aa^j  exp  (a* a) |0) 
=  exp  (aaf)  |0)  exp  (-|a|2/2) 


(3.27) 
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where  the  Baker-Cambell-Hausdorff  theorem  is  used  to  get  the  second  line  and  the 
property  of  the  annihilation  operator  on  the  lowest  weight  state  is  used  to  get  to  the 
third  line. 

The  normalization  exp  |a|2/2^  can  be  transferred  to  the  metric  defining  in  this 
way  a  coherent  state  that  is  not  normalized,  but  which  has  the  advantage  of  being  analytic 
in  a.  The  canonical,  analytic  coherent  state  is: 

\a)  =  exp  (aaf)  |0)    .  (3.28) 


In  general  I  will  work  with  analytic  coherent  states.  This  choice  provides  powerful 
mathematical  properties,  as  can  be  noted  from  the  fact  that  the  three  most  common 
time-dependent  variational  principles  are  equivalent  [72]  with  such  trial  wavefunctions. 

A  general  multiconfigurational  wavefunction  for  the  electrons  has  been  formulated 
by  Weiner,  Deumens  and  Ohrn  [83]  using  the  mathematics  of  vector  coherent  states  [84, 
85].  However,  I  will  not  describe  that  formulation  here  but  restrict  the  details  to  the  case 
of  the  single  determinantal  approximation  for  the  electronic  wavefunction,  which  is  the 
model  that  has  been  developed  into  a  working  program. 

The  coherent  state  description  of  the  single  determinant  leads  to  the  Thouless 
parametrization  [73].  Assuming  an  orthonormal  spin-orbital  basis  set  made  up  of  K 
vectors  and  N  particles  to  be  described  by  a  single  determinant,  this  determinant  can  be 
found  as  a  unitary  transformation  acting  on  the  orbital  expansion  coefficients  of  some 
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reference  single  determinant  |^o)-  I  will  define  this  reference  state  as 

N 

i*o)=n6/iuac>  <3-29> 
/=i 

where  b]  (6,)  are  the  fermion  creation  (annihilation)  operators  of  the  basis  set  used.  The 
reference  state  is  a  lowest  weight  state  for  the  irreducible  representation  [1N0(A:  ^]  of  the 
group  U(K)  with  generators  b]bj  because  the  weight  operators  6-6,  have  eigenvalue  1  for 
1=1... JV,  and  0  for  i=N+l,..X.  The  stability  group  of  the  reference  state  is  U (N)  ®U(K- 
N).  Introducing  a  complex  parametrization  of  the  coset  space  U (K)/U(N)  ®  U(K  —  N) 
an  element  g  of  the  group  U(k)  is  written  as  g  =  ch  as  previously  discussed.  Extending 
U(K)  to  GL(K,  C)  and  using  the  Gauss  factorization  of  c  (see  notational  remarks  below) 


this  leads  to  a  parametrization  of  the  coset  space  in  terms  of  the  complex  parameters  z, 
with  x  and  y  to  be  determined  as  functions  of  z  from  the  condition  that  g  is  unitary. 

In  the  last  expression  and  in  what  follows  I  use  a  notation  intended  to  make  the 
equations  simpler  to  follow.  In  the  atomic  basis,  there  is  a  "hole"  subspace  generated  by 
the  first  N  orbitals  which  make  up  the  reference  single  determinant  and  the  "particle" 
space  made  up  of  the  rest  of  the  basis  spin-orbitals.  The  field  operators  of  the  basis 
states  which  belong  to  the  holes  will  be  written  with  a  bullet,  i.e.  b'\  b',  and  those 
associated  to  the  particle  states  will  have  an  open  circle  b„,  6°.  In  the  same  way, 
when  referring  to  the  molecular  orbitals  made  up  from  the  atomic  orbitals,  the  occupied 
molecular  orbitals  have  a  bullet  and  the  virtual  orbitals  have  open  circles.  So  bullets  will 
always  be  used  to  describe  hole  states  in  a  given  basis  and  open  circles  will  represent  the 


(3.30) 
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particle  states.  Matrices  associated  with  the  subspaces  of  holes  or  particles  will  also  carry 
these  symbols.  Matrix  blocks  which  carry  a  prime  or  a  double  prime  indicate  the  upper  or 
lower  off-diagonal  block  respectively.  Matrices  associated  with  the  complete  one-particle 
space  will  not  carry  any  symbols.  So,  for  example,  we  may  write  the  identity  matrix 
as  /  =  (7Q  7°0) 

Let  T  be  the  unitary  irreducible  representation  of  U(K)  in  Fermi-Fock  space.  The 
coherent  state  is  defined  as 

|*,)  =T(c  fc)|*o) 

-T(c)T(fc)|»0) 

=r(<OI*o> 

=aT((72#    7°0))|*o>  (3-31) 

=aU\^  +  E  6;S«W> 

i=l  \         j=N+l  J 

/  N      K  \ 
=aexp £  E  !••>• 

\i=l  j=JV+l  / 

The  second  step  uses  the  fact  that  h  is  in  the  stability  group  of  the  reference  state,  and 
that  the  rightmost  factor  of  c  in  (3.30)  modifies  only  the  virtual  space  and  leaves  the 
reference  state  unaffected.  The  middle  factor  of  c  only  gives  a  constant  a  when  acting 
on  the  reference  state.  The  constant  a  can  be  determined  from  the  normalization  of 
the  coherent  state.  Note  that  the  representation  is  considered  as  acting  on  the  orbital 
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coefficients,  not  on  the  basis  functions.  This  is  in  accordance  with  the  active  point  of 
view  of  coordinate  transformations.  I  work  with  the  unnormalized  coherent  state 


This  is  the  Thouless  representation  of  a  determinantal  wave  function  with  the  elements 
of  the  (K  —  N)  x  N  matrix  z  as  time-dependent  parameters. 

Thouless  [73]  worked  in  another  direction  to  get  to  this  result.  Because  his  derivation 
clarifies  the  meaning  of  the  z  parameters  I  will  include  it  here.  If  U  is  the  matrix 
which  transforms  from  the  basis  spin-orbitals  to  the  occupied  and  virtual  orbitals  of  the 
fermionic  system,  then 


A' 


(3.32) 


t=i 


(3.33) 
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Following  Thouless, 

m=l[c?\vac) 
t=i 

1=1  \j=l  j=N+l 

=  ft  (  E  K  +  E  E W'  I  <«  )  W  (3-34) 
,=i  \/=i  \        ;=#+!  *=i 

=°n(V+  E  DjVsr'M 

t=l  \  j=N+lk=l 

«=1  \       jWV+1  jb=l  /  /=1 

where  the  invariance,  up  to  a  constant  a,  of  a  determinantal  wave  function  under  a  linear 

transformation  of  its  occupied  spin  orbitals  is  used.  Then  the  coherent  state  is  recovered 

as  follows 

i*>=n(l+  e  *;vni*o) 

i=l  \       j=N+l  ) 
N  K 


=n  n  (i+*-v?)i*o) 

i=l  >=JV+1 
JV  A" 

=n  n  exp^-^ivo) 

t=l  j=N+l 

=exp(E  E      1 1*>, 

\i=l  ]=N+1 


(3.35) 
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where  X)i=i  UflPhi    =  zji-      Eq-  (3.35)  we  have  used  the  nilpotency  of  operators 

The  determinantal  wavef unction  is  expressed  as  det{x(*/)}  where 

K 

Xi  =  1>i+  £  (1  <  i  <  AT).  (3.36) 

are  nonorthogonal  (but  linearly  independent)  dynamical  orbitals.   The  corresponding 

unoccupied  dynamical  orbitals  may  be  chosen  as 

N 

X;  =  ^-E2;.^     (N+l<j<K)  (3.37) 
t=i 

and  although  mutually  nonorthogonal  they  are  orthogonal  to  the  occupied  space. 
The  occupied  dynamical  orbitals  can  then  be  expressed  as 

(P    *°)(rz)  (3-38) 

i.e.  a  partitioned  row  vector  of  basis  spin  orbitals  times  a  partitioned  rectangular  matrix. 
Consider 


-z    1°  )(  z      1°  )  "I  -z  +  z 


-z*  +  z*\ 
1°  +  zz*  J 


_(P  +  zU  0 
~  V      0         1°  +  zz^ 

This  shows  that 


(3.39) 


(3.40) 


defines  dynamic  occupied  and  virtual  orbitals  that  span  orthogonal  spaces.  It  follows  that 
the  projectors  on  these  spaces  add  up  to  the  unit  operator. 

/  -  pocc  =  (3.41) 
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These  relations  are  used  many  times  when  simplifying  various  expressions  in  the  evolution 
equations. 

All  of  the  above  assumed  a  basis  of  orthonormal  spin-orbitals.  In  practical  applica- 
tions it  is  preferable  to  work  in  the  nonorthogonal  atomic  spin-orbital  basis.  In  this  case, 
since  the  anticommutation  relations  of  the  field  operators  is  not  the  canonical  one  but 
includes  the  overlap  matrix  of  the  orbitals,  it  is  not  possible  to  write  the  coherent  state 
as  an  exponential.  However,  it  is  still  possible  to  express  the  molecular  orbitals  in  the 
same  way  as  in  Eq.  (3.36)  with  modified  z's: 


The  orbitals  of  the  virtual  space  are  not  as  simple  to  express  as  in  the  orthonormal 
basis.  Assuming  a  form 


K 


Xi  -  *  +  E  (!<»<*). 


(3.42) 


j=N+l 


(3.43) 


:=1 


such  that  the  (K-N)xN  matrix  v  satisfies 


A*  +  ztA't  +  a'z  +  z\  A°z  +  A'  +  z+A'V  +  z*A° 

uA*  +  vA'z  +  A'*  +  A°z      A0  +  vA'  +  A'V  +  uA,ut 


) 


(3.44) 
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i.e.  the  two  subspaces  are  orthogonal,  then 

vA*  +  vA'z  =  -A,f  -  A°z  (3.45) 

and  hence 

v  =  -  (A't  +  A°z)(A,  +  A'z)"1 

(3.46) 

„t  =  _  (a*  +  2f A't) _1  (A'  +  *t A0). 

The  Equations  of  Motion 

The  preceding  tools  will  now  be  applied  to  find  specific  equations.  These  equations 
are  the  ones  that  are  written  into  the  computer  code  DYNAMO.  Here  I  will  derive  the 
general  form  of  the  equations  and  comment  on  their  properties.  The  details  can  be  found 
in  the  appendices. 

Three  approximations  are  made: 

1.  In  the  TDVP  the  limit  of  narrow  Gaussian  wavepackets  is  taken.   The  nuclear 
parameters  are  then  the  position  and  momentum  of  the  nuclei. 

2.  The  electrons  are  described  by  means  of  a  single  determinant. 

3.  The  spin-orbital  electronic  basis  is  truncated  to  a  finite  number  of  orbitals  without 
electron  translation  factors. 

The  truncated  basis  limitation  is  always  imposed  if  an  actual  calculation  is  to  be 
performed.  The  lack  of  ETF's  is  something  which  will  be  remedied  in  the  near  future 
to  allow  us  to  analyze  collisions  at  high  speeds.  For  lower  speeds  associated  with  most 


41 

chemical  reactions,  the  velocities  involved  are  sufficiently  small  and  ETF's  are  not  as 
critical.  The  single  determinantal  wavefunction  can  be  improved  upon  by  including  more 
configurations.  The  theory  to  do  this  for  the  electronic  states  can  be  found  in  Ref.  83. 
For  nuclei,  there  is  the  question  of  finding  a  good  quantum  representation  for  them. 
Presently  work  is  being  done  in  this  direction.  A  first  approach  is  to  describe  them  via 
frozen  Gaussians  without  taking  the  limit  of  narrow  wavepackets.  These  are  the  canonical 
coherent  states  discussed  in  the  previous  section.  However,  in  work  done  with  Deumens 
and  Ohrn  [3],  we  have  shown  that  a  description  in  which  the  nuclear  wavefunctions 
are  not  given  the  freedom  to  split,  the  dynamics  becomes  nearly  identical  to  that  done 
using  classical  nuclei.  Currently,  work  is  in  progress  to  describe  nuclei  going  beyond 
the  Gaussian  wavepacket  approximation. 

Initially  I  will  assume  that  the  electronic  basis  is  complete  and  is  independent  of 
the  positions  of  the  nuclei.  Then  I  will  show  how  to  derive  equations  for  two  different 
electronic  basis  sets  which  are  more  adequate  for  calculations.  The  new  equations  can 
be  obtained  via  a  symplectic  transformation. 

The  basis  set  that  does  not  depend  on  the  nuclear  positions  will  be  referred  to  as 
a  Non-Following  Basis  (NFB).  Another  basis  of  interest  is  one  with  orbitals  that  are 
instantaneously  centered  at  the  position  of  the  nuclei,  but  which  do  not  have  associated 
with  them  the  motion  of  the  nuclei,  i.e.  they  do  not  contain  ETF's.  This  basis  set  will 
be  referred  to  as  the  Static-Following  Basis  (SFB).  A  basis  set  centered  at  the  position  of 
the  nuclei  with  ETF's  will  be  called  the  Dynamic  Following  Basis  (DFB).  For  practical 
calculations  with  a  truncated  basis  it  is  desirable  to  use  either  the  SFB  or  the  DFB. 
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Otherwise,  after  some  time,  the  nuclei  may  be  far  from  the  vicinity  of  the  basis  set 
location,  which  would  be  inadequate  for  a  correct  description  of  the  system. 

There  are  two  equivalent  ways  to  derive  the  TDVP  equations  for  narrow  Gaussians. 
One  is  to  find  the  equations  of  motion  and  then  take  the  limit  of  zero  width  and  the 
other  is  to  take  the  limit  before  actually  varying  the  action.  I  will  do  the  latter.  In  this 
case  the  action  becomes 


A-/[iE(A-*-A-*Hi: 


Mi  +  (z\z) 


dt  (3.47) 


— * 

where  H  is  the  molecular  Hamiltonian,  including  the  nuclear  repulsion  term,  Rj  and 
Pj  are  the  canonical  positions  and  momenta  of  the  nuclei  and  \z)  is  the  coherent  state 


description  of  the  electronic  single  determinant  as  defined  in  Eq.  (3.32). 


Defining 


Zk  =  Rk+  iPk 


(3.48) 


we  can  write 


1 


(3.49) 
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The  electronic  overlap  matrix  is  given  by 

S(z'\z)=(z'\z) 

N  N 
«=1  ;=1 

=  det((W+  £  £  (3.50) 

=  det(^tJ  +   £  «fa-4>)v) 

V  ifc=./V+l  / 

=  det(/'  +  z'tz) 

and  if  in  the  derivation  of  the  TDVP  equations  the  quantity  5  is  replaced  by 


S(('\()Sd(Z'*,Z)  (3.50) 

then  Eq  (3.11)  can  be  used  to  get  the  equations  of  motion  for  the  electronic  and  nuclear 
parameters: 


dE 

dz* 


dE 

(3.51) 
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dE 
dza 
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dZk 
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where  the  unit  matrix  comes  from  Sc\.  In  matrix  form,  these  equations  become 
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which,  with  the  real  form  of  the  nuclear  coordinates,  can  be  expressed  as 


\ 

f  z\ 

(dE/dz*\ 

z* 

dE/dz 

it 

dE/dR 

/ 

\  dE/dP  J 

(3.52) 


(3.53) 


A  transformation  that  leaves  a  set  of  dynamical  equations  invariant,  i.e.  does  not 
change  the  metric  at  all,  is  called  a  symplectic  transformation  and  if  the  metric  contains 
only  ones  and  zeroes,  the  transformation  is  said  to  be  canonical.  We  are  interested  in  a 
generalized  symplectic  transformation  that  leaves  the  structure  of  the  equations  invariant, 
but  changes  the  values  of  the  matrix  elements  of  the  metric.  The  invariant  is  the  Poisson 
bracket.  First  we  consider  the  transformation  from  a  NFB  to  a  DFB.  It  has  the  general 

form  . 

z  =  z{z,  Z,  Z  ) 

(3.54) 

Z  =  Z{Z) 

or  with  the  real  form  of  the  nuclear  coordinates 

z  =  z{z,R,  P) 


R  =  R(R) 


(3.55) 


P  =  P{P). 

With  this  notation  I  indicate  that  the  z  parameters  are  independent  dynamical  variables, 
i.e.  they  are  independent  of  the  "new"  nuclear  coordinates  and  momenta.  In  order  to 
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transform  the  Poisson  bracket,  we  need  the  matrix  of  partial  derivatives,  Jacobian,  J 

(3.56) 
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such  that 


J  = 


I '  c*  0  0  0\ 

0  c  0  0 

r*  r  I  0 

\p*  P  0  // 


(3.57) 


where  c  =  dz/dz,  r  =  dz/dR,  and  p  =  dz/dP.  Using  the  fact  that  the  inverse  of  a 
Jacobian  is  the  Jacobian  of  the  inverse  transformation  on  the  Poisson  bracket  results  in 
the  transformation 

{fyg}=difM-1dg 

=  {/,<?}- 

=  d]fJ]M~lJdg 


(3.58) 


of  the  total  phase  space  metric  M  =  JMJ\  with 

{ic*CcT        0  ic*CrT 

0  -icC*J  -tcCM 

ir*CcT  -irC*c*  ir*CrT  -  irC*r* 

\ ip*CcT  -iPC*c*  I  +  ip*CrT  -  ipC*r* 


M  = 


ic*CpT  \ 
-icC*pi 
-I  +  ir*CpT  -  irC'pi 
ip*CpT  -  ipC*p*  ) 


(3.59) 


We  then  obtain  for  the  dynamical  equations  in  the  new  coordinates  the  expression 
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Here  we  define  the  matrices 


Cxk  = 


d2\nS(z*,R!,P',z,R,P) 
dz*dXu 


R'=R,P'=P 


and 


d2\nS(z*,R',P',~z,R,P) 
dXLdYi 


R'=R,P'=P 


(3.61) 


(3.62) 


CXkY,  =  -21m- 
with  X  and  Y  standing  for  R  or  P. 

For  the  case  of  coherent  states  in  terms  of  the  atomic  basis  centered  on  the  nuclei  the 
calculation  of  the  metric  involves  the  overlap  of  two  coherent  states  with  different  nuclear 
geometries.  If,  furthermore,  electron  translation  factors  are  included  in  the  orbitals,  the 
overlap  also  depends  on  the  velocity,  and  hence  the  momentum,  of  the  nuclei.  We  find 

S(z'*,     P\  z,  R,  P)  =  det  (A*  +  A'z  +  *'f  A"  +  z'^A°z)  (3.63) 
where  the  overlap  matrix 

A(R!,P',R,  P)  (3.64) 
depends  on  the  nuclear  positions  and  momenta.  We  have  dropped  the  tilde  on  the  z 
parameters,  something  we  do  from  now  on,  since  we  will  only  deal  with  a  following 
basis.  Note  that  the  overlap  matrix  of  two  different  bases  is  not  Hermitian,  and  that  it 
becomes  the  unit  matrix  when  the  two  nuclear  configurations  and  momenta  coincide. 

In  the  general  derivation  it  is  assumed  that  the  basis  set  can  depend  on  both  the  nuclear 
positions  and  momenta.  This  is  the  case  for  the  DFB.  If  the  dependence  of  the  basis  on 
the  nuclear  parameters  is  only  through  the  nuclear  positions,  the  equations  become 
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These  are  the  equations  that  have  been  coded  into  DYNAMO,  using  the  results  found 
in  Appendix  B  for  a  nonorthogonal  basis 

The  interpretation  of  these  equations  is  very  interesting.  I  will  indicate  in  a  general 
form  the  most  salient  features.  The  first  point  to  make  is  that  the  equations  are  time- 
reversible.  This  comes  from  the  time-reversibility  of  the  TDVP. 

Energy  is  conserved.  This  can  be  seen  most  easily  using  the  Poisson  brackets 

E  =  {E,E}  =  0  (3.66) 

Another  property  that  can  be  checked  directly  is  conservation  of  total  momentum  of 
the  system.  The  expectation  value  of  the  electronic  momentum  is  pe\  =  (z\pei\z)  /  (z\z) 
where  the  operator  in  brackets  is  the  total  electronic  momentum  operator.  The  time 
derivative  of  this  quantity  is 

ft,  =  £  (2ImC^5  -  CW4  -  CrpPi)  (3.67) 

since  the  time  dependence  only  comes  from  the  electronic  and  nuclear  parameters.  By 
inspecting  Eq.  (3.60),  multiplying  the  third  row  of  the  metric  and  adding  over  all  the 
nuclei  gives 

£  {-2lmCnj + cnaA  +  c*p'p)  -Y.'pi  =  Y.jwE  <3-68> 
*=i  k=\    k=\  °ni 

Because  the  basis  set  is  centered  at  the  position  of  the  nuclei,  the  energy  is  invariant  to 
a  global  translation  of  the  nuclei  and  so  the  right  hand  side  vanishes.  Using  the  last  two 
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expressions  we  get: 


Pel  +  ^Tot  =  0 


(3.69) 


which  shows  the  conservation  of  total  momentum. 

For  a  truncated  NFB  the  preceding  result  is  not  valid  since  in  that  case  it  can  be 
shown  that 


but  in  a  limited  basis,  Ehrenfest's  theorem  is  not  generally  satisfied,  which  implies  that 


Thus,  in  a  truncated  non-following  basis  total  momentum  would  not  be  conserved.  This 
property  is  only  satisfied  in  the  basis  sets  which  exhibit  following. 

The  next  question  I  consider  is  the  interpretation  of  the  electronic  equations  of 
motion.  The  meaning  is  identical  when  the  nonorthogonal  atomic  basis  is  used,  but  the 
more  obscure  mathematics  coming  from  the  nonorthogonality  makes  interpretation  more 
difficult.  Using  the  orthogonal  basis  and  drawing  from  the  results  found  in  Appendix  A 
to  write  down  the  equation  of  motion  for  the  z's 


(3.70) 
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where  the  vertical  bar  after  the  overlap  matrix  indicates  that  we  are  taking  the  derivative 
with  respect  to  the  unprimed  coordinates  in  Eq.  (3.64).  This  gives  the  following 
expression  for  the  forces  in  matrix  form 

iz  +  7°)(AvAA|  +  ftVAA| 

(3.73) 

which  can  be  further  reduced  to 

H  +  i^i-z    /^(^V^AI  +  ^V^AI 

=  (-»   /»)((.  + TV(V.Mir) J  (J2*)  (3-74) 

-(-« ^(c)- 

The  right  hand  side  has  a  very  simple  interpretation:  the  Fock  matrix  acts  on  the  occupied 
states  and  is  then  projected  on  the  virtual  space  and  only  a  nonzero  projection  on  the 
virtual  states  gives  rise  to  a  change  in  the  z  coefficients.  The  Fock  matrix  in  a  basis 
which  contains  the  ETF's  of  the  form  exp(iv  •  r)  contains  two  extra  terms  in  the  h  part 
of  the  Fock  matrix  derived  from  the  operation  of  the  electronic  kinetic  energy  operator 
on  the  ETF's.  One  of  them  is  associated  with  the  kinetic  energy  v2/2,  and  the  other 
exactly  cancels  the  term  coupling  with  the  nuclear  velocity  in  the  metric.  Thus,  only 
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a  term  associated  with  the  nuclear  acceleration  is  left  If  we  assume,  for  a  moment, 
that  we  are  treating  only  a  single  atom  and  that  the  electron  is  initially  in  its  ground 
state  and  that  for  some  reason  the  nucleus  is  being  accelerated,  then  the  Fock  matrix 
will  not  be  generating  any  coupling,  but  the  term  with  the  nuclear  acceleration  generates 
a  dipole  coupling  between  the  basis  functions  which  couples  to  the  acceleration  of  the 
nucleus.  This  generates  dynamics  in  the  electrons.  What  happens  is  that  the  electrons 
are  being  excited  through  the  nuclear  acceleration.  These  are  the  couplings  I  mention  in 
the  previous  chapter  which  I  have  not  seen  implemented  in  the  literature. 

If  a  basis  without  ETF's  is  used,  then  the  Fock  matrix  reverts  to  the  standard  time- 
independent  form  and  the  coupling  to  the  nuclear  velocity  in  the  metric  does  not  get 
cancelled.  Analyzing  this  term  in  a  complete  basis  and  assuming  that  the  nucleus  moves 
with  a  constant  velocity  shows  that  this  term  in  the  metric  is  associated  with  the  translation 
operator.  After  a  time  St  the  nucleus  has  moved  SR  =  R  -  St.  The  electronic  states  are 
to  be  described  in  terms  of  a  basis  centered  at  the  new  position.  The  old  basis  functions 
must  be  expressed  in  terms  of  the  new  ones,  so  the  transformation  matrix: 


should  be  evaluated  to  apply  on  the  expansion  coefficients.  For  a  small  translation  the 
electron  translation  operator  is 


and  its  matrix  elements  in  a  basis  which  is  centered  on  the  nucleus  (so  that  derivatives 


r  _  SR)\<j>(r- Rj)  =  (fl exp  (sii  ■  Vr-) \<j>) 


(3.75) 


(3.76) 
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with  respect  to  the  electron  and  nuclear  coordinates  can  be  interchanged)  gives 


6ij  +  ifiStV  AA\ 


(3.77) 


Because  we  use  an  unnormalized  state,  the  changes  in  the  z's  are  computed  as  projected 
from  the  occupied  to  the  unoccupied  states  through  the  translation  operator.  Dividing  by 
St  leads  to  a  time  rate  of  change  of  the  z's  which  is  precisely  the  second  term  found 
in  Eq  3.74.  So  the  action  of  this  coupling  is  to  take  care  of  the  fact  that  the  electron 
coefficients  are  being  continually  defined  with  respect  to  a  basis  which  has  been  translated. 
This  coupling,  which  is  sometimes  called  "unphysical"  [6],  is  just  ensuring  the  correct 
description  of  the  electronic  state  in  a  translating  basis. 


Tools  must  be  developed  to  understand  the  results  of  a  time-dependent  calculation. 
I  have  developed  two  such  tools  which  will  be  described  here.  The  first  one  is  used  to 
analyze  state-to-state  transfer  probabilities  for  the  p+H  collision.  The  results  found  in 
the  static  following  basis  must  be  projected  on  eigenstates  of  the  moving  hydrogen  atom. 
The  eigenstates,  which  were  described  in  Chapter  2  when  discussing  ETF's,  are 


Interpretive  Tools 


Am(f-  R)  =  exp  (iv-r-  i(v2/2  +  tnlm)t)4>nim(r  -  R^j 


(3.78) 


where  D  indicates  the  dynamic  nature  of  the  orbital.  The  transformation  from  the  static 
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to  the  dynamic  basis  is  determined  from: 

i*>  =  E  i^>a* 
k 

(3.79) 

=  £  l^lja! 
/.it 

=  £  itf  >a? 
/.*r 

It  is  easily  seen  that  the  part  in  the  dynamic  phase  which  is  multiplied  by  the  time 
t  is  cancelled  in  the  projection  Yl\(f>F)(<f>F\  smce  *s  independent  of  the  electronic 
coordinates.  What  is  required,  therefore,  is  the  transformation  matrix  defined  by 

(exp(iv-  f)(f>nim\<f>n>vm>)  =  {<j>nlm\exp(-iv  •  r)\<f>n,i>ml)  (3.80) 

which  when  applied  to  the  expansion  coefficients  in  the  SFB  gives  the  dynamic  expansion 
coefficients.  These  will  be  called  the  boosted  coefficients. 

To  compute  such  a  transformation  is  not  trivial  in  a  hydrogenic  or  STO  basis. 
However,  if  the  basis  set  is  written  in  terms  of  Gaussians,  the  computation  can  be  done 
analytically.  DYNAMO  relies  on  the  integral  package  ABACUS  written  by  Helgaker, 
Jensen  and  Jorgensen  [86]  using  Gaussian  basis  functions.  The  basis  functions  used  in  the 
calculations  of  the  p+H  collision  involve  the  hydrogen  Is,  2s  and  2p  functions  written 
as  contractions  of  Gaussians. 

The  s-  and  p-type  normalized  Gaussian  functions  have  the  form 

<f>3(f;a)  =  (2a/7r)3/4  exp  (— ar2) 

(3.81) 

<f>px,{f;a)  =  x,(l28a5/7r3)1/4exp  (-ar2) 
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Then,  the  fundamental  matrix  elements  from  which  to  construct  the  boost  transformation 
matrix  are 


3 /2 

W*(")l«tp(-«?-»W«OT) 
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vxvy    (2(a(3f2\5/2       (       v2  \ 
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(3.82) 

and  from  the  fixed  linear  combinations  of  Gaussians  describing  the  atomic  basis  it  is  a 
simple  matter  to  evaluate  the  projection  from  the  SFB  to  the  boosted  basis. 

These  results  are  implemented  in  a  program  called  BOOST.  It  uses  the  output  from 
DYNAMO,  performs  the  transformations  and  computes  state-to-state  probabilities. 
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For  a  many-electron  molecular  process  in  which  the  final  product  is  made  up  of 
several  distinct  fragments  we  want  to  analyze  the  output  of  a  DYNAMO  run  and  find 
the  probability  of  having  a  given  number  of  electrons  in  each  fragment. 

In  general  we  will  also  want  state-to-state  probabilities,  even  for  the  case  of  a 
single  molecule,  not  many  fragments.  For  a  many-electron  system  these  states  are  time- 
independent  ones  which  approximate  the  eigenstates  of  the  fragments  or  molecule.  I  have 
developed  the  following  method  to  obtain  such  information  from  the  single  determinantal 
output  from  DYNAMO. 

As  discussed  in  the  previous  section,  our  starting  point  is  a  single  determinant  made 
up  of  molecular  spin-orbitals  (MSO's)  which  are  a  linear  combination  of  nonorthogonal 
atomic  spin-orbitals  (ASO's).  This  determinant  is  some  final  or  even  intermediate  result 
of  a  DYNAMO  run.  I  will  call  the  MO's  that  are  obtained  in  this  way  dynamic  MSO's. 
The  determinant  is 


We  can  interpret  the  final  result  in  terms  of  time-independent  multiconfigurational 
states  (I  will  use  MC  to  denote  multiconfigurational,  not  Monte  Carlo)  which  are  an 
approximation  to  the  exact  eigenstates  of  the  system.  If  there  are  several  fragments 
then  we  will  want  to  interpret  the  result  in  terms  of  MC  states  corresponding  to  each 
moiety.  The  simplest  example  is  that  of  a  collision  process:  a  fragment  hits  another, 
a  reaction  occurs  and  then  fragments  separate.  Another  might  be  two  systems  coming 


A' 


(3.83) 


j=N+l 
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in  close  contact,  some  reaction  occurring  and  then  separating.  In  all  these  cases  there 
are  two  aspects: 

1.  There  are  different  possible  reaction  channels,  i.e.  a  different  number  of  electrons 
associated  with  each  fragment. 

2.  We  want  to  understand  the  chemistry  in  terms  of  "local"  states,  where  by  "local"  I 
mean  molecular  states  written  only  in  terms  of  the  ASO's  of  the  fragment  of  interest. 

The  general  procedure  is  to  do  a  MC  calculation  for  each  of  the  fragments  and  for  all 
the  channels  of  interest  within  a  fragment.  This  will  become  clearer  in  a  moment.  For  a 
given  fragment  and  channel  (where  by  channel  I  mean  a  given  number  of  electrons  in  the 
fragment)  a  standard  SCF  is  performed  to  determine  a  set  of  MSO's  which  are  used  to 
perform  a  MC  calculation  to  a  given  level  of  interest.  These  I  will  call  reference  MSO's 
for  the  fragment  and  channel.  In  principle,  it  is  not  necessary  to  use  a  HF  state  as  the 
reference  state.  A  more  general  reference  state  can  be  used  if  needed.  One  possibility 
is  to  use  natural  orbitals. 

What  is  needed  is  the  transformation  between  the  reference  MSO's  and  the  nonorthog- 
onal  basis  functions.  Calling  this  transformation  U,  we  have 


<t>j  =  XiUij 


(3.84) 


where  <f>  are  the  reference  MSO's,  and  the  relation 


USUI  / 


(3.85) 


is  used,  with  S  the  overlap  matrix  of  the  ASO's. 
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Next  we  need  to  expand  the  state  we  have  in  terms  of  the  fragments.  To  make  this 
clearer  assume  there  are  only  two  fragments.  It  is  straightforward  to  generalize  to  more 
fragments.  The  fragments  contain  a  well  defined  subset  of  basis  ASO's.  The  ASO's 
in  each  fragment  ideally  have  no  overlap  with  those  of  the  other  fragment.  If  there  is 
some  overlap  I  assume  it  is  small  and  can  be  ignored  for  the  calculation.  Each  dynamic 
MSO  is  a  sum  of  the  form 

W  =  S  x>  c>«  +  ^  Xl  Cli  =  ^  +  ^  (3'86) 

jcA  UB 

The  single  determinant  is  then  written,  using  the  fact  that  a  determinant  is  multilinear  in 
its  columns  or  rows,  as  (we  also  normalize  the  total  wavef unction  to  1  here) 

=  [{\<Pi<f4-v»1t\}  +  {ivf  vi  ...v?j}l 

(3.87) 

+  +  •••}  + ...  +  {lA?-4l}]/v»> 

where  each  brace  contains  {j^j  single  determinants  corresponding  to  all  the  possible 
ways  of  getting  (iV  -  M)  distinct  states  of  fragment  A  and  M  distinct  states  of  fragment 
B,  (M  =  0, 1...A0.  Another  way  of  understanding  these  states  is  that  they  are  all 
the  possible  determinants  with  (./V  —  M)  electrons  in  fragment  A  and  M  electrons  in 
fragment  B.  In  general,  most  of  these  will  be  zero,  or  very  unimportant  compared  to 
others.  For  example,  if  a  fragment  with  L  electrons  comes  in  contact  with  another  of 
M  electrons,  there  might  be  one  electron  transferred  from  A  to  B,  so  that  the  relevant 
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terms  of  the  expansion  are  only  those  with  either  L  electrons  in  A  or  L  -  1  electrons 
in  A.  Otherwise,  there  will  be  2^  determinants  in  the  preceding  expansion,  which  might 
be  necessary  in  some  cases  (such  as  in  the  a  +  He  collision  in  which  one  determinant 
represents  2-electron  transfer,  another  represents  2-electron  excitation  and  the  other  two 
determinants  represent  only  1 -electron  transfer). 

The  scheme  just  described  sets  the  stage  for  the  MC  calculations.  According  to  the 
channels  of  interest  in  the  system,  a  MC  calculation  is  performed  to  find  the  ground  and 
excited  states  for  the  fragments.  Perhaps  at  most  a  CI  Singles  and  Doubles  (CISD)  will 
be  possible  due  to  computational  limitations  in  the  size  of  the  system,  and  even  then  it 
might  be  truncated  to  only  a  certain  number  of  excited  states.  This  choice  will  depend 
on  the  system  being  studied. 

The  next  step  is  to  expand  the  determinants  of  Eq.  (3.87)  in  terms  of  the  MSO 
transformation  U  for  each  fragment  and  appropriate  channel.  Redefining  U  as  Ufo, 
where  C  refers  to  the  fragment  and  M  to  the  number  of  electrons  in  the  fragment,  and 
defining  the  reference  MSO's  of  fragment  C  with  M-electrons  as  <^.t,  the  dynamic 
MSO's  are  written  as 


where  Eq.  (3.84)  is  used  to  write  the  dynamic  MSO  in  terms  of  the  reference  MSO's. 
Any  determinant  of  Eq.  (3.87  )  can  be  written  in  terms  of  determinants  of  the  reference 


(3.88) 
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MSO's,  for  example 

il,S2,  — ,tjv 

These  determinants  are  in  terms  of  basis  functions  that  are  the  same  as  those  used  in  the 
MC  calculations.  Since  the  MC  calculation  is  done  for  the  fragments,  the  states  we  want 
to  compare  with  are  made  up  with  those  MC  states  in  the  following  manner: 

Ol  .•••>£} 

where  the  fragment  MC  states  are 

(3.91) 

*X;j     =    £    l^i;ir-^f;iil*Oi,..iL}i  • 
0i>  -it} 

Since  I  assume  that  the  ASO's  in  A  and  B  are  orthogonal,  the  overlap  between  the 
reference  MSO's  of  different  fragments  are  null.  When  an  overlap  is  taken  between  these 
MC  states  and  the  dynamic  determinant  expanded  in  terms  of  the  reference  MSO's,  what 
is  obtained  is  a  sum  of  determinants  of  overlap  matrices.  Because  of  the  property  of  zero 
overlap  between  MSO's  of  different  fragments,  the  overlap  determinants  can  always  be 
brought  into  a  block  form,  which  is  just  the  product  of  the  determinants  of  each  fragment. 

In  this  way  we  come  to  the  final  result.  The  overlap  determinants  are  either  (-l)p  (to 
account  for  any  orbital  ordering  within  the  determinants)  or  0  due  to  the  orthonormality 
of  the  MSO's. 
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With  this  approach,  questions  which  can  be  asked  include: 

1.  Which  MC  state  is  the  most  important  in  the  final  state. 

2.  What  is  the  probability  of  finding  the  system  in  a  given  MC  state. 

Some  probability  will  always  be  lost  to  higher  MC  states  that  are  not  taken  into  account, 
unless  a  full  CI  is  performed  and  all  possible  channels  retained.  By  using  common  sense 
and  physical  and  chemical  intuition,  the  size  of  this  loss  should  be  kept  reasonably  small. 
It  is  also  easy  to  check  how  much  probability  is  lost  to  eigenstates  not  considered  by 
adding  the  calculated  probabilities  and  subtracting  the  sum  from  the  total  of  one. 

This  procedure  has  been  coded  in  a  program  called  PROJECT.  It  does  not  incorporate 
the  state-to-state  analysis  through  a  MC  projection  yet.  What  it  does  do  is  calculate  the 
probability  of  ending  up  in  different  channels.  It  is  used  in  the  work  presented  in  the 
next  chapter  to  calculate  the  1-  and  2-electron  transfer  probabilities  in  a  +  He. 


CHAPTER  4 
RESULTS 

The  p+H  Collision 

Introduction 

Proton-hydrogen  collisions  have  been  the  subject  of  abundant  experimental  studies  in 
the  last  20  years.  Thus  they  provide  an  excellent  test  for  new  time-dependent  theoretical 
approaches  to  electron-nuclear  dynamics. 

Most  theoretical  work  on  this  system  is  performed  at  collision  energies  above  1  keV  in 
the  lab  frame.  Calculations  use  either  straight  line  or  Coulomb  trajectories  for  the  nuclei. 
Such  approaches  yield  electron  transfer  and  excitation  cross  sections  in  agreement  with 
experimental  results,  even  if  strict  conservation  laws  of  energy  and  total  momentum  are 
violated.  However,  for  lower  energies  the  motion  of  the  protons  is  sufficiently  slow  that 
different  trajectories  significantly  alter  the  results.  It  is  then  necessary  to  treat  the  electron- 
nuclear  and  nuclear-nuclear  interactions  correctly  throughout  the  collision  process. 

Runge  et  al  [60]  use  the  eikonal  approximation  to  treat  the  nuclei  and  solve  the  time- 
dependent  electronic  density  matrices  in  the  linearized  TDHF  approximation  as  described 
in  Chapter  2.  They  analyze  proton-hydrogen  collisions  in  the  10  eV-1  keV  collisional 
energy  range  using  only  Is  states  in  the  basis.  They  also  investigate  the  effect  of  using 
straight  line  trajectories,  Coulomb  trajectories  and  effective  potentials.  They  show  the 
importance  of  using  the  correct  trajectory  for  energies  below  1  keV. 

60 
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Fritsch  and  Lin  [53]  as  well  as  Kimura  and  Lane  [56]  have  recently  published  reviews 
on  the  semiclassical  close-coupling  method,  and  examined  the  proton-hydrogen  collision 
in  some  detail,  as  well  as  some  other  collisions.  This  method,  as  well  as  some  other 
methods  they  mention,  all  use  prescribed  trajectories. 

Fritsch  and  Lin  use  the  close  coupling  scheme  with  an  extended  basis  which  they  call 
the  AO+  method  [54,  55].  It  consists  of  a  basis  of  Is,  2s,  2p  (and  in  Ref.  55,  also  n=3 
states)  atomic  orbitals  corresponding  to  the  free  hydrogen  atom  as  well  as  atomic  orbitals 
corresponding  to  the  united  atom  (He).  The  He  orbitals  are  expected  to  be  important  for 
small  impact  parameters.  Their  work  explores  the  1-75  keV  energy  range. 

The  numerical  integration  of  the  time-dependent  Schrbdinger  equation  is  realized  by 
Grim  et  al  [87]  at  2  keV  using  a  numerical  integration  on  a  grid,  and  avoiding  the  use 
of  basis  sets. 

Other  work  includes  that  by  Liidde  and  Dreizler  [88,  89]  who  solve  the  time  dependent 
Schrodinger  equation  using  a  large  set  of  Hylleraas  type  basis  functions. 

Few  theoretical  studies  are  done  for  energies  below  1  keV.  The  most  recent  by 
Hunter  and  Kuriyan  [49,  50]  for  collision  energies  between  .0001  eV  and  10  eV,  uses 
the  PSS  method  to  separate  the  nuclear  from  electronic  degrees  of  freedom.  Davis  and 
Thorson  published  a  similar  study,  but  for  an  energy  range  going  from  0  to  0.2  eV  [51] 
and  correcting  for  some  errors  in  the  Hunter  and  Kuriyan  work.  The  results  for  the 
two  studies  in  the  overlapping  energy  range  are  very  similar.  I  will  be  comparing  with 
Hunter  and  Kuriyan's  work  since  the  information  they  report  is  in  the  form  of  tables, 
while  Davis  and  Thorson  only  publish  graphs  with  a  low  precision.  In  both  cases  the 
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nuclei  are  treated  quantum  mechanically,  and  only  the  molecular  \sa  and  2pa  states 
are  retained.  Twenty  years  before  this  work,  Dalgarno  and  Yadav  [66]  also  used  a  the 
PSS  approach  to  treat  this  problem  for  energies  starting  as  low  as  0.25  eV  and  going 
up  to  energies  of  100  keV.  A  little  before  that  Bates  and  Dalgarno  studied  this  problem 
using  the  Born  approximation,  for  the  energy  range  0-250  keV  [90].  The  results  of  this 
early  work  shows  that  a  perturbative  scheme  like  the  Born  approximation  works  well  for 
energies  above  10  keV,  but  fails  for  lower  energies  by  several  orders  of  magnitude. 

An  important  development  at  energies  above  1  keV  is  the  inclusion  of  electron 
translation  factors  (ETF's)  in  the  basis.  In  the  limit  of  a  complete  set,  such  factors  are 
not  needed,  but  the  question  always  lingers  as  to  how  large  a  basis  must  be  to  account 
correctly  for  the  couplings.  A  review  of  these  factors  can  be  found  in  Chapter  2.  When 
using  a  molecular  orbital  basis,  different  ETF's  have  led  to  different  results,  particularly 
for  the  25  excitation  and  transfer  cross  section.  This  can  be  seen,  e.g.  in  the  10  basis 
close-coupling  calculations  of  Crothers  and  Hughes  [91]  and  Kimura  and  Thorson  [52]. 

I  use  the  computer  code  DYNAMO  to  calculate  several  properties  of  the  proton- 
hydrogen  collision  as  a  first  test  of  the  general  method.  No  ETF's  are  used,  which  limits 
the  collisional  energies  that  can  be  investigated.  I  choose  an  energy  range  between 
0.02-4000  eV  which  spans  the  very  low  energy  regime  investigated  by  Hunter  and 
Kuriyan  and  Davis  and  Thorson  [49-51],  as  well  as  the  higher  energies  of  other  research 
groups.  The  transfer  or  excitation  probabilities  to  different  orbitals  are  calculated  by 
first  projecting  the  results  on  a  moving  hydrogenic  basis  (i.e.  one  with  the  factor  etmv  T 
included),  as  explained  in  the  previous  chapter. 
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The  momentum  conservation  law  discussed  in  Chapter  3  allows  DYNAMO  to  be 
used  to  show  explicitly  the  interchange  of  momentum  between  the  three  particles  of  this 
system. 
Results 

Approximate  hydrogenic  Is,  2s,  and  2p  states  are  expressed  in  terms  of  6  Gaussians 
each.  The  Is  and  2p  states  are  taken  from  Stewart  [92].  The  hydrogen  2s  was  optimized 
by  me.  The  coefficients  and  exponents  used  can  be  found  in  Table  4-1. 


Table  4-1:  Contraction  coefficients  c  and  exponents  a  for  the  basis  used  in  this  work. 


Is  orbital* 

2s  orbital 

2p  orbital8 

c 

c 

a 

c 

a 

1.30334 

6.51095 

-3.75318 

1.21392 

1.01708 

1.214225 

xlO1 

xlO"2 

xlO1 

xlO"2 

xlO1 

xlO"2 

4.16492 

1.58088 

-7.61767 

2.67784 

4.25860 

2.64900 

xlO1 

xlO1 

xlO1 

xlO"2 

xlO1 

xlO2 

3.70563 

4.07099 

3.22382 

1.60950 

4.18036 

6.099425 

xlO1 

xlO1 

xlO"2 

xlO1 

xlO1 

xlO"2 

1.68538 

1.18506 

1.77665 

4.71398 

1.73897 

1.585355 

xlO1 

xlO1 

xlO1 

xlO1 

xlO1 

4.23592 

4.93615 

5.42463 

1.68965 

3.76794 

5.10090 

xlO"2 

xlO2 

xlO"2 

xlO1 

9.16360 

2.31030 

1.01853 

9.22099 

3.75970 

2.577175 

xlO"3 

xlO1 

xlO2 

xlO"3 

a:  Ref.  [92]. 


The  hydrogen  atom  in  its  ground  state  is  initially  placed  at  rest  at  the  origin  of  the 
laboratory  coordinate  system  and  the  proton  at  a  distance  of  50  a.u.  from  the  origin  and 
given  different  impact  parameters  and  an  initial  momentum  corresponding  to  the  energy 
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of  the  collision.  The  system  is  allowed  to  evolve  until  the  nuclei  are  again  separated 
by  50  a.u. 

For  collision  energies  from  10  eV  to  4  keV,  the  approximate  n=l  and  n=2  hydrogenic 
basis  is  used.  A  comparison  between  this  larger  basis  and  a  smaller  basis  of  a  single  Is 
state  per  center  is  made  at  10  eV  and  the  difference  in  total  cross-section  is  found  to  be 
less  than  2%.  Thus  for  energies  from  0.02  eV  to  10  eV  only  the  Is  functions  are  used. 
Impact  parameters  are  chosen  differently  for  various  energy  ranges:  from  0.1  to  7.9  a.u 
in  steps  of  0.2  a.u.'s  for  energies  from  100  eV  to  4  keV;  from  0.1  to  9.9  a.u.'s  in  steps 
of  0.2  aoi.'s  for  energies  of  10  eV  to  40  eV,  and  0.0  to  13.9  a.u.'s  in  steps  of  0.1  a.u.'s 
for  energies  of  0.02  eV  and  0.18  eV.  When  only  the  Is  basis  is  used,  the  initial  and  final 
separations  were  taken  to  be  30  a.u.'s. 

The  transition  probabilities  are  computed  after  projecting  on  the  basis  with  ETF's 
included.  Then  transfer  and  excitation  cross- sections  are  computed  from: 


where  n,l  define  the  state,  and  P(E,  b)  is  the  probability  at  the  given  collision  energy 
E  and  impact  parameter,  b. 

As  a  first  step  to  determine  if  DYNAMO  is  capable  of  reproducing  good  results 
for  this  system  I  plot  in  Fig.  4-1  the  probability  of  electron  transfer  times  the  impact 
parameter  as  a  function  of  the  impact  parameter  and  compare  it  to  an  exact  integration 
of  the  time-dependent  Schrodinger  equation  where  the  wavefunction  is  represented  on  a 
grid.  This  last  work  is  done  by  Griin  et  al  [87]  at  2  keV.  Also  shown  are  results  of  work 
done  by  Liidde  and  Dreizler  who  use  a  close  coupling  method  without  ETF's  and  a  basis 


(4.1) 
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of  64  Hylleraas-type  basis  functions  [88].  The  agreement  between  our  results  and  those 
of  Griin  et  al  is  excellent,  while  Ludde  and  Dreizler's  results  are  not  that  good.  Work  by 
Fritsch  and  Lin  using  the  close-coupling  with  ETF's  and  their  AO+  method  also  shows 
very  good  agreement  with  the  results  by  Griin  et  al  [53]. 
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Figure  4-1:  Weighted  transition  probabilities  for  total  electron  transfer  at  2  keV  as  a  function  of 
impact  parameter.  All  data  in  atomic  units.  Dots  are  the  results  by  Griin  et  al.  The  full  line 
are  the  results  of  calculations  by  Ludde  and  Dreizler.  The  dash  dotted  line  are  the  results  from 
DYNAMO.  (See  text  for  references). 


Table  4-2  lists  the  total  transfer  cross-sections  for  the  p+H  collisions  for  the  energies 
studied  and  compares  them  with  the  experimental  results. 

Fig  4-2  plots  the  experimental  total  transfer  cross-sections  and  our  calculations.  From 
the  previous  table  and  this  plot,  it  is  clear  that  our  results  reproduce  experiments  over 
five  orders  of  magnitude  in  the  kinetic  energy  of  the  colliding  proton. 
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Table  4-2:  Total  transfer  cross-sections  for  proton  colliding  with  a  hydrogen  atom  (x  10~16  cm2). 


Collision  Energy  (eV).        Total  Transfer  Experiment  (Energy,  eV) 

Cross-section. 


0.02 

69.55 

U.  10 

j  / 

59.4      .9/7(0.  lo)^ 

in  n 

^7 

37.0  3  _4.7(9.7)a'a 

20 

33.22 

32.8±4.8  (22.2)^ 

40 

30.25 

30.0±4.9  (40.4)a-d 

100 

25.60 

23.7±3.5  (109.6)8 

500 

19.44 

18.9±3.2  (500)b 

1000 

16.78 

16.3±2.9  (1000)b 

2000 

14.07 

13.9±3.5  (2000)b 

3000 

12.43 

12.1+0.61  (3040)c 

4000 

11.33 

11.1±0.55  (3820)c 

a:  Newman  et  ai,  Ref.  [93]. 
b:  Gealy  and  Van  Zyl,  Ref.  [94]. 
c:  McClure,  Ref.  [95]. 

d:  At  energies  below  100  ev,  collisions  are  between  protons  and  deuteron  atoms. 


Newman  et  al.  [93]  measure  total  transfer  cross-sections  for  energies  between  0.18 
eV  and  300eV.  For  energies  below  100  eV,  they  use  protons  colliding  on  deuterium 
atoms.  Hunter  and  Kuriyan  [49,  50]  use  a  partial  wave  expansion  for  the  nuclei  moving 
on  a  Born-Oppenheimer  lsa,  2pa  electronic  potential  energy  surfaces  at  low  energies 
(below  10  eV).  Their  results  indicate  that  below  1  eV  appreciable  differences  appear 
between  the  total  transfer  cross-sections  of  a  proton  on  deuterium  and  that  of  the  proton 
on  hydrogen.  As  the  former  process  is  slightly  endothermic  ( AE  =  0.0037  eV),  this  result 
is  to  be  expected  at  some  point.  We  study  collisions  between  protons  and  deuterium  and 
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D+-hydrogen  collisions.  Our  results,  as  well  as  Hunter  and  Kuriyan's,  for  energies  below 
1  eV  are  shown  in  Table  4-3.  The  calculations  with  the  deuterium  atom  and  ion  are  done 
using  the  same  basis  as  the  proton  on  hydrogen,  except  that  the  exponents  in  the  basis 
are  changed  to  reflect  the  change  in  the  reduced  mass  of  the  deuteron-electron  system. 


Energy  (eV) 

Full  circles:  Ref.  93.  Full  squares:  Ref  94.  Full  triangles:  Ref.  95. 

Figure  4-2:  Total  transfer  cross  sections  from  0.02  eV  to  4000  eV.  Comparison  of  experiment  to 
theory.  DYNAMO  computations  are  the  larger  open  circles  joined  by  solid  line. 

Comparing  the  results  from  DYNAMO  with  those  of  Hunter  and  Kuriyan  [49,  50] 
and  noting  that  the  experimental  value  of  the  cross-section  at  0.18  ±  0.06  eV  for  the 
p+D->H+D+  is  59.4  IJ^6  x  10"16  cm2  [93],  their  result,  at  0.2  eV,  is  just  within  the  error 
bar  and  considerably  lower  than  the  one  reported  here.  The  present  calculations  show 
only  a  slight  difference  between  the  proton  colliding  against  a  hydrogen  atom  and  against 
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a  deuterium  atom.  Even  at  0.02  eV  the  difference  is  not  significant.  Since  at  0.18  eV 
DYNAMO  agrees  more  closely  with  the  experimental  value,  it  would  be  interesting  if 
experiments  could  be  repeated  at  these  and  perhaps  even  lower  energies  to  verify  when 
and  if  the  isotopic  effect  becomes  important.  Due  to  obvious  problems  in  colliding  a 
proton  to  a  hydrogen  atom,  it  would  be  easier  and  more  interesting  to  measure  collisions 
of  a  proton  on  deuterium  and  collisions  of  D+  on  hydrogen.  Hunter  and  Kuriyan  report 
significant  differences  in  the  transfer  cross-sections  of  these  two  processes  for  energies 
below  0.2  eV. 


Table  4-3:  Total  transfer  cross-sections  for  collisions  of  a  proton  on  hydrogen  and  deuterium 
atoms,  as  well  as  that  of  a  D+  ion  on  hydrogen,  for  energies  below  1  eV,  compared  to  results  by 
Hunter  and  Kuriyan  (H&K)  (Cross-section  units  xlO  _16  cm  2). 


Energy  (eV) 

p+H^H+p 

p+D-+H+D+ 

D++H->D+p 

H&Ka 

0.02 

72.63 

60.78 

72.17 

Present  work 

0.02 

69.55 

69.23 

69.25 

H&Ka 

0.2 

55.88 

50.64 

51.31 

Present  work 

0.18 

57.32 

56.89 

a:  Refs.  49,  50. 


For  transfer  and  excitation  cross-sections  to  n=2  states  at  energies  of  10-40  eV,  the 
z  coefficients  to  the  n=l  states  are  below  the  accuracy  requested  from  the  integrator. 
Thus  transfer  and  excitation  cross- sections  to  these  states  are  reported  here  only  for 
energies  of  100  eV  and  above.  The  results  are  shown  in  Table  4-4.  Fig.  4-3  and  Fig. 
4-4  compare  total  2p  and  25  cross-sections  with  experiment  for  excitation  and  transfer 
processes  respectively. 


69 

Table  4-4:  Excitation  and  transfer  cross-sections  for  2s,  2p\,  2pz  states  and  total  2p  cross-sections 
(x  1016  cm2  ). 


Enerev 

Excitation 

Transfer 

2s 

2px 

2pz 

2p 

2s 

2px 

2pz 

2p 

lOOeV 

0.01 

0.42 

0.72 

1.14 

0.19 

0.56 

0.68 

1.24 

500eV 

0.98 

12.45 

2.64 

15.09 

2.16 

11.35 

2.57 

13.92 

lOOOeV 

1.59 

23.76 

3.54 

27.30 

4.61 

11.28 

9.38 

20.66 

2000eV 

6.78 

30.42 

3.71 

34.13 

10.37 

10.65 

18.98 

29.63 

3000eV 

10.54 

31.77 

1.85 

33.62 

12.99 

10.78 

32.66 

43.44 

4000eV 

11.51 

33.13 

1.90 

35.02 

13.64 

10.48 

45.22 

55.70 

Fig.  4-3  and  Fig.  4-4  compare  total  2p  and  2s  cross-sections  with  experiment  for 
excitation  and  transfer  processes  respectively. 

The  results  from  DYNAMO  are  very  close  to  experimental  values  for  total  2p 
excitations.  They  are  slightly  larger  than  experiments  after  2  keV.  These  results  and 
other  recent  theoretical  work  by  others  coincide  with  the  more  recent  measurements 
rather  than  with  the  older  Stebbings  et  al  [96]  results,  which,  for  excitations,  appear  to 
be  too  low  and  for  transfer  seem  to  be  too  high. 

The  calculated  2s  excitation  cross-sections  fall  right  between  the  values  of  the  two 
experiments  performed  at  these  energies.  Other  theoretical  treatments  [54,  55,  89,  52] 
locate  the  2s  excitation  cross-sections  below  these  experimental  values. 

The  lack  of  ETF's  for  excitations  is  not  as  important  as  in  the  case  of  electron  transfer 
since  little  momentum  is  transferred  to  the  target  for  the  important  impact  parameters. 
The  only  problem  is  that  the  projectile  states  are  not  well  reproduced,  since  they  lack 
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important  ETF's,  an  omission  which  limits  the  accuracy  of  the  excitation  cross- sections 
for  higher  energies. 


4.00 


0.00 


0.00  1.00  2.00  3.00 

Energy  (keV) 


4.00 


5.00 


Open  diamonds:  Ref.  97.  Open  boxes:  Ref.  98.  Open  triangles:  Ref.  96.  Solid  boxes:  Ref.  99. 
Solid  triangles:  Ref.  100. 

Figure  4-3:  n=2  excitation  cross  sections.  Comparison  of  DYNAMO  results  with  experiments. 
Total  2p  excitation  calculations  are  shown  by  open  circles  joined  by  a  solid  line.  2s  excitation 
calculations  are  shown  by  solid  circles  joined  by  a  dotted  line. 

DYNAMO  2p  transfer  cross-sections  agree  with  experiments  up  to  about  2keV. 
Above  that  they  are  too  large.  For  2s  transfer  the  results  are  well  above  the  existing 
experimental  values.  Other  theoretical  studies,  in  particular  by  Kimura  and  Thorson  [52] 
and  by  Fritsch  and  Lin  [54,  55]  do  a  better  job  of  reproducing  experimental  results.  Our 
omission  of  ETF's  in  these  calculations  is  probably  the  cause  of  the  discrepancies  at 
these  collisional  energies. 
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Figure  4-4:  n-2  transfer  cross  sections.  The  symbols  and  lines  have  the  same  meaning  as  in 
the  previous  figure. 

The  integral  alignment  provides  a  measure  of  the  relative  excitation  to  the  different 
2p  states  of  the  target  atom.  The  available  experimental  results  are  from  Hippler  and 
collaborators  [101,  102].  The  measurements  they  perform  start  from  1  keV  and  go  up 
to  energies  of  25  keV. 

The  integral  alignment  is  computed  as: 


A2o  = 


<Tq  +  2(7! 


(4.2) 


where,  for  our  calculations: 


<70  =  crPz 


(4.3) 


2<T1  =  °px 

The  z-axis  is  the  initial  direction  of  motion  of  the  proton. 
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Fig.  4-5  compares  my  results  to  experiment.  In  Refs.  101,  102  the  experiments  are 
compared  to  theoretical  work  by  Liidde  and  Dreizler  [89]  and  Fritsch  and  Lin  [54,  55].  As 
the  energy  drops  to  1  keV,  all  these  theoretical  results  converge  to  an  alignment  of  46%, 
just  above  the  experimental  error  bars,  and  none  predict  the  decrease  in  the  alignment 
below  2  keV  that  experiments  show.  The  straight  line  approximation  used  in  these  other 
works  appears  to  cause  errors  in  the  alignment  at  these  lower  energies  which  are  not 
present  when  using  a  full  electronic-nuclear  dynamics  treatment  as  in  DYNAMO.  My 
results  reproduce  the  decrease  in  the  alignment  and  predict  a  sharp  fall  in  the  alignment 
as  energy  continues  to  drop.  This  warrants  further  experimental  work. 


xlO1 


0.00        1.00        2.00        3.00        4.00        5.00  6.00 

Energy  (keV) 


Solid  triangles:  Ref.  101.  Solid  boxes:  Ref  102. 

Figure  4-5:  Computed  integral  alignment  and  experimental  results,  as  a  function  of  collision 
energy. 
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Finally,  Fig.  4-6  shows  the  x  components  of  the  nuclear  and  electronic  momenta 
which  are  initially  zero  for  a  collision  at  40  eV  and  an  impact  parameter  of  3.0  a.u.  The 
basis  set  used  for  this  calculation  is  a  Is  orbital  on  each  nucleus.  The  x  component  of 
the  total  momentum  remains  zero  throughout  the  calculation. 


Momentum  (a.u.) 

4.0 

3.0 

2.0 

1.0 

0.0 
•1J0 
-2.0 
•3.0 
-4.0 

350.         400.         450.         500.         550.         600.  650. 
Tim*  (a.u.) 

Figure  4-6:  Expectation  value  of  x  component  of  electron  momentum  and  nuclear  momentum 
as  a  function  of  time.  Middle  curve  is  the  electron  expectation  value,  the  upper  curve  is  the 
momentum  corresponding  to  the  incoming  proton  and  the  lower  curve  is  that  of  the  proton 
originally  at  rest  at  the  origin. 
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The  a+He  Collision 


Introduction 

After  finding  the  excellent  agreement  with  experiment  that  DYNAMO  can  achieve 
for  a  one-electron  system,  it  is  important  to  determine  how  it  can  handle  a  system  with 
more  than  one  electron.  One  possibility  is  to  test  it  on  the  p+He  collision.  However  this 
system  is  essentially  a  one-electron  transfer  process  and  it  seems  more  important  to  test 
our  method  on  a  system  for  which  different  channels  are  possible.  The  simplest  system 
for  which  there  is  enough  experimental  work  is  the  a+He  collision. 

Experiments  to  determine  the  one-  and  two-electron  transfer  on  the  a+He  collisions 
have  been  performed  by  several  groups  [103-108].  The  most  comprehensive  work  has 
been  done  by  Afrosimov  et  al  in  the  middle  70's  [103,  104].  They  used  a  coincidence 
measurement  technique  to  draw  information  about  different  excited  states  in  the  products. 
They  used  3  He  in  their  experiments.  The  results  by  this  group  and  results  by  Bayfield 
and  Khayrallah  [105]  and  by  Berkner  et  al  [106]  all  agree  within  the  error  bars  of  about 
20%  for  each  of  these  groups.  Except  for  the  results  by  Afrosimov  et  al  and  Berkner 
et  al,  all  of  the  experiments  are  done  at  collision  energies  above  10  keV.  Berkner  et  al 
provide  results  at  energies  above  7  keV,  and  Afrosimov  et  al  present  results  at  energies 
above  2.5  keV.  Experiments  at  lower  energy  have  been  performed  by  Hertel  and  Koski 
[108]  for  energies  above  0.5  keV  and  by  Latypov  et  al  [107]  for  energies  above  0.1 
keV.  These  latter  results  have  such  large  uncertainties  that  is  is  fruidess  to  compare 
calculations  with  them. 
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On  the  theoretical  side,  work  on  this  system  has  been  done  using  the  PSS  approach 
[109],  different  close-coupling  methods  [110,  111]  and  using  the  TDHF  [112,  113].  An 
investigation  of  variational  improvement  of  the  TDHF  is  carried  out  by  Gazdy  and  Micha 
[114-1 17]  for  the  a+He  collision.  However,  this  work  is  formal,  carried  out  for  the  one- 
electron  transfer  at  energies  between  30-100  keV  and  using  only  a  minimal  basis  set 
of  Is  orbitals. 

The  PSS  work  covered  an  energy  range  of  10-100  keV  [109]  obtaining  close 
agreement  with  experiment.  The  basis  used  is  an  expansion  in  7  single  determinants 
each  determinant  formed  with  one-electron  diatomic  orbitals  (i.e.  orbitals  found  for  the 
diatom  containing  only  a  single  electron)  and  not  using  ETF's.  The  close-coupling  work 
by  Kimura  [110]  explores  an  energy  range  from  2-400  keV  using  15  determinants  from 
STO  basis  functions.  He  uses  an  MO  approach  with  MO-ETF's.  His  results  also  agree 
well  with  experiment,  using  straight  line  or  Coulomb  trajectories  depending  on  the  energy. 
Gramlich  ex  al  [111]  use  a  Gaussian  orbital  expansion  with  ETF's  to  look  at  the  energy 
range  of  8-^400  keV.  They  only  use  straight  line  trajectories.  This  affects  the  quality  of 
their  results  at  lower  energies.  They  use  a  total  of  29  determinants  for  their  calculations. 
At  higher  energies,  where  curved  trajectories  are  not  as  important,  their  results  are  closer 
to  experiments  than  those  of  Kimura. 

The  TDHF  work  is  done  in  the  energy  range  of  30-150  keV  by  Devi  and  Garcia  [112] 
and  20-160  keV  by  Stich  et  al  [113].  The  first  used  a  prescribed  Coulomb  trajectory  and 
several  simplifying  assumptions  on  the  TDHF  equations.  In  the  second  case  straight  lines 
were  prescribed  for  the  nuclear  motion.  The  electronic  basis  functions  were  expanded 
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on  a  large  basis  of  Hylleraas  type  functions  (  the  actual  number  varied  according  to 
the  distance  between  nuclei,  but  was  on  the  order  of  150).  The  results  of  these  two 
calculations  are  good.  However  for  energies  around  30  keV  and  under  the  one-electron 
transfer  is  overestimated  in  both  cases.  This  has  been  explained  as  a  lack  of  electronic 
correlation,  i.e.  the  use  of  a  single  determinant.  It  has  also  been  stated  by  Stich  et  al  that 
to  lift  the  single  determinantal  description  is  a  very  difficult  problem. 

Here  I  use  the  single  determinant  approximation  for  the  electrons  at  a  much  lower 
energy  range  than  the  TDHF  results  cited  above.  The  purpose  was  to  determine  the 
importance  of  different  basis  sets  as  well  as  to  test  a  single  determinantal  description 
at  lower  energies.  I  look  at  the  energy  range  of  4-10  keV  and  find  the  total  one-  and 
two-electron  cross-sections. 
Results 

Four  different  basis  sets  are  used.  The  first  one  (Basis  I)  is  identical  to  the  one  used 
for  the  p+H  collision,  except  that  the  exponential  coefficients  are  multiplied  by  a  factor 
of  3.9967  to  account  for  the  charge  of  +2  for  the  He  nucleus  and  also  for  the  effect  of  the 
mass  of  the  He  atom  on  the  Bohr  radius.  The  second  basis  (Basis  II)  is  made  up  of  two 
contracted  s  functions,  the  first  made  up  of  3  Gaussians  and  the  second  of  1  Gaussian  and 
a  p  set  made  from  a  single  contracted  Gaussian  [118].  This  basis  is  a  basis  optimized  for 
hydrogen.  The  third  basis  (Basis  III)  is  similar  to  Basis  II,  but  optimized  for  helium.  The 
fourth  basis  (Basis  IV)  has  4  s  basis  functions  from  the  contraction  of  7  s-type  Gaussians 
and  2  p  basis  functions  made  from  the  contraction  of  3  p-type  Gaussians  [119].  This  last 
basis  set  is  also  optimized  for  the  He  atom. 
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The  calculations  are  started  with  the  nuclei  separated  by  50  a.u.  with  the  He  atom 
at  the  origin  of  coordinates  and  in  the  ground  state  as  determined  by  the  optimization  of 
the  He  atom  in  each  basis.  The  ion  is  given  an  impact  parameter  ranging  from  0.0  to  2.9 
a.u.  in  steps  of  0.1  a.u.  for  collision  energies  of  2.5,  4,  6,  8  and  10  keV,  although  all 
these  energies  are  not  used  for  all  basis  sets.  At  the  end  of  a  run  the  program  PROJECT, 
described  in  Chapter  3,  is  used  to  determine  the  probability  for  one-  and  two-electron 
transfer.  The  transfer  cross  section  is  calculated  using 


for  A  being  one  of  the  two  transfer  processes. 

The  range  of  0.0-2.9  a.u.  is  used  because  at  2.9  a.u.  the  one-electron  transfer 
probability  is  of  the  order  of  10"5,  while  for  the  2-electron  transfer  it  is  10"10. 

The  results  of  the  calculations  are  tabulated  in  Table  4-5  for  one-electron  transfer 
and  Table  4-6  for  two-electron  transfer.  Experimental  results  of  Afrosimov  [103,  104] 
are  also  shown. 

Table  4-5:  One-electron  transfer  cross  sections  using  three  different  basis  sets  (see  text  for 
description  of  the  basis  sets)  are  compared  with  experiments  by  Afrosimov.  (Cross-section 
units  x  10"17  cm2.) 

Energy        Basis  I  Basis  II         Basis  UJ        Basis  IV        Exp.  (error) 


(4.4) 


2.5  keV 


4.54 


3.49 


4.0  (0.6) 
4.0  (0.6) 
3.9  (0.6) 
3.9  (0.6) 
3.9  (0.6) 


4  kev 


4.97 


6.32 


6.77 


6keV 


7.56 


7.58 


7.76 


8.36 


8keV 


14.99 


lOkeV 


10.60 


10.34 
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Table  4-6:  Two-electron  transfer  cross  sections  using  three  different  basis  sets  (see  text  for 
description  of  the  basis  sets)  are  compared  with  experiments  by  Afrosimov.  (Cross-section 
units  x  10"17  cm2.) 


Energy 

Basis  I 

Basis  II 

Basis  III 

Basis  IV 

Exp.  (error) 

2.5  keV 

18.4 

15.3 

31.0  (4.7) 

4keV 

15.0 

16.1 

14.8 

25.0  (3.7) 

6keV 

15.2 

17.2 

16.9 

15.6 

23.0  (3.5) 

8keV 

15.0 

21.0  (3.1) 

lOkeV 

16.0 

17.9 

20.0  (3.0) 

The  results  for  the  different  basis  sets  used  are  very  similar,  with  slightly  better  results 
for  the  two-electron  transfer  using  the  basis  set  Basis  IL  The  nearly  constant  experimental 
cross- section  for  the  two-electron  transfer  is  qualitatively  reproduced  by  DYNAMO. 

The  one-electron  transfer  is  not  reproduced  very  well.  At  the  low  end  (2.5  keV) 
the  agreement  is  very  close.  Otherwise,  the  cross-sections  are  off  by  a  factor  of  2  or 
3.  This  agrees  with  overestimation  of  one-electron  transfer  probabilities  found  using  the 
TDHF  by  Devi  and  Garcia  [112]  and  Stich  et  al  [113],  although  they  never  got  as  low 
as  this  in  energy. 

The  two-electron  transfer  process  is  essentially  a  resonant  transfer  phenomenon  [103]. 
Electrons  tend  to  transfer  from  ground  state  to  ground  state  of  the  He  atom.  On  the  other 
hand,  the  one-electron  process  tends  to  be  a  transfer  to  an  excited  n=2  or  higher  state  of 
He+  [104].  By  using  a  single  determinant,  we  force  it  to  describe  both  these  processes 
in  a  single  configuration.  This  could  be  a  possible  explanation  for  the  results  not  being 
completely  satisfactory  at  these  energies. 
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As  the  energy  goes  down,  the  single  determinant  seems  to  provide  a  better  picture 
for  the  one-electron  transfer.  This  could  be  indicating  another  possible  reason  for  a  poor 
description  of  one-electron  transfer:  the  lack  of  ETF's  in  the  basis  is  skewing  couplings 
and  degrading  the  results. 

The  method  yields  acceptable  results.  Being  a  factor  of  2  or  3  off  from  experiment 
is  not  bad!  The  results  have  the  right  orders  of  magnitude  and  certain  trends  are  repro- 
duced. Having  tested  four  different  basis  sets  shows  that  the  results  have  approximately 
converged  within  the  approximations  of  the  method.  The  differences  with  experiment  do 
not  come  from  too  small  a  basis  set. 

These  results  as  well  as  the  ones  for  the  p+H  collision  indicate  that  ETF's  should  be 
included  if  DYNAMO  is  to  be  used  to  predict  results  of  collisions  in  the  keV  range. 

Another  conclusion  is  that  the  single  determinant  may  not  be  able  to  predict  quan- 
titatively accurate  results  when  singlet  and  triplet  states  both  play  significant  roles  in  a 
chemical  process,  such  as  for  the  a+He  collision.  This  requires  an  extension,  which  in 
our  case,  contrary  to  the  problems  faced  by  Stich  et  al,  is  quite  simple  to  implement. 


CHAPTER  5 
CONCLUSIONS 


The  single  determinantal  Hartree-Fock  approximation  in  electronic  structure  calcula- 
tions is  not  always  very  accurate.  Yet  it  is  usually  an  excellent  starting  place  for  methods 
which  include  electronic  correlation  and  it  is  quite  capable  of  generating  qualitatively  and 
even  quantitative  results  which  are  in  good  agreement  with  experiment  [81].  For  large 
systems  the  Hartree-Fock  approach  is  the  starting  point  for  even  further  simplifications 
leading  to  semiempirical  methods. 

In  the  same  way,  a  single  determinantal  approach  to  time-dependent  methods  seems 
to  be  a  good  starting  point.  Although  this  assertion  needs  to  be  investigated  further,  it 
is  clear  from  the  results  presented  here  and  those  we  have  already  published  [4]  that  the 
single  determinantal  approach  does  yield  quite  good  results. 

The  alternatives  at  present  are  to  either  do  dynamics  on  potential  energy  surfaces  using 
quantum  mechanical  methods,  or  using  a  full  CI  approach  to  time  dependent  methods. 
The  first  approach  can  produce  excellent  results  for  three  atom  systems  when  only  two 
or  three  potential  energy  surfaces  are  important  in  the  dynamics.  Not  many  chemical 
systems  fall  in  this  category,  however.  The  full  CI  approach  is  embodied  in  the  close- 
coupling  methods.  These  methods,  as  implemented  today,  ignore  nuclear  dynamics  and 
so  are  limited  to  collisions  in  the  keV  ranges.  By  doing  this  they  violate  conservation 
of  energy,  momentum  and  angular  momentum.  They  also  face  the  same  problem  that 
CI  does  in  electronic  structure  calculations:  computer  time  and  memory  needed  grow 
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factorially  with  the  number  of  electrons  and  basis  functions.  For  this  reason  close- 
coupling  calculations  have  been  limited  to  systems  of  no  more  than  2  electrons.  An 
intermediate  position,  i.e.  one  where  a  systematic  multiconfigurational  approach  is  used 
without  going  to  full  CI  has  not  been  formulated  in  the  literature. 

The  field  of  time-dependent  dynamics  with  a  single  determinant  is  not  new  in  physics, 
although  its  application  to  chemical  systems  is  fairly  recent.  The  advantage  of  the  END 
equations  presented  here  as  compared  to  similar  methods  lies  in  the  transparency  of 
its  approximations  and  the  ease  with  which  they  can  be  improved  upon.  This  method 
has  already  been  formulated  in  such  a  way  as  to  use  a  multiconfigurational  electronic 
description  [83],  although  a  computer  code  has  not  been  implemented  to  work  with  such 
a  state  yet.  The  conservation  properties  of  the  equations  are  also  important.  No  other 
single  determinantal  method  has  been  developed  to  the  extent  of  the  one  presented  here. 

For  the  two  systems  studied  here  I  have  shown  that  the  lack  of  ETF's  can  hinder 
accurate  results  for  energies  above  a  few  keV.  If  collisions  or  reactions  are  to  be  studied 
at  lower  energies  than  this  then  the  method  yields  good  results.  We  also  cover  a  range 
of  energies  that  methods  based  on  PES  or  the  close-coupling  approach  can  handle,  i.e. 
between  10  eV  and  1  keV. 

The  single  determinant  approach  must  be  tested  further  to  learn  more  about  it's 
strengths  and  weaknesses.  A  first  upgrade  is  to  write  a  new  integral  code  which  is  more 
adequate  for  the  kind  of  operations  that  a  time-dependent  method  requires.  For  example 
it  should  include  electron  translation  factors,  and  it  should  be  vectorizable  to  make  it 
efficient  on  modern  computers. 
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Another  extension  is  in  two  very  different  directions;  on  one  hand  to  include  more 
determinantal  states  for  chemical  processes  where  electron  correlation  is  important,  and 
on  the  other  hand  to  develop  a  semiempirical  approximation  starting  from  a  single 
determinant  to  permit  the  study  of  large  systems. 

Another  line  of  investigation  is  to  include  a  quantum  description  of  the  nuclei  in  such 
a  way  as  to  permit  wavepacket  splitting  for  the  nuclei.  The  approach  to  take  is  one  in 
which  not  all  nuclei  have  a  quantum  description  since  heavier  nuclei  are  less  likely  to 
experience  significant  quantum  effects  from  tunneling  or  zero  point  energy  fluctuations. 


APPENDIX  A 

DETAILS  OF  DERIVATION  IN  ORTHONORMAL  BASIS 

The  derivations  in  this  appendix  and  the  next  can  be  found  in  the  technical  report 
"Coherent  State  Approach  to  Time  Evolution  with  a  Hartree-Fock  State"  by  Deumens, 
Diz  and  Ohm  [120]. 

In  Chapter  3  the  overlap  kernel  for  the  electrons  is  found  to  be: 

S{z'\  z)  =  (z'\z)  =  det(7"  +  z'fz)  (A.l) 

from  which  we  can  derive  expressions  for  its  derivatives 

d 


dzi 


-ln(zV)|^=2=(.|z)-1^^MvMinor(M),^|2^ 


(A.2) 


((/•  +  A)-V 


jk 

In  the  first  step  the  determinant  is  expanded  in  the  j'-th  row  with  M  denoting  the  matrix 

M  =  r  +  z*z.  (A.3) 

The  next  step  uses  the  expression  for  the  inverse  of  a  matrix  in  terms  of  its  minors 

Minor(M)- 

v       >J%        det(M)  V  ' 

83 


84 

and  the  fact  that  none  of  the  elements  z/y  do  occur  in  the  minor  ij  because  it  does  not 
contain  column  /  of  M. 

Similarly  we  find  for  the  derivative  with  respect  to  the  complex  conjugate  parameter 


Hi 


(A.5) 


-(.(r+A)-") 


This  can  be  written  in  matrix  form  as 

9Tln(z'|Z)|z-=z=  (/'  + A)"1^ 


dzT 

^ln(z'|z)|2.=*  =  *(/'  + A)"*. 


(A.6) 


dz 

Eq.  (3.63)  shows 

S{z'\B!,P\z,R,P)  =  (z'\B!,P'\z,R,P) 

(A.7) 

=  det(A*  +  A'z  +  z'^A"  +  z'^A°z) 
from  which  the  derivative  with  respect  to  the  nuclear  position  can  be  evaluated.  It  is 

the  purely  imaginary  quantity 

V*  \n(z',R',P'\z,R,P)  I 

R"    x    '     '     1  1  \z'=z,R'=R,P'=P 

(A.8) 

=  Tr(/-  +  ^^(v^A-l  +  VRA'\z  +  z^&A"\  +  z^nA°\z), 

where  we  used  the  following  property 

dlndet  (M)      -    .  ddet  (M)  dMtJ 

dx       =d6t(M)  ?"lj^T" 

«i 


=  TrM" 


5t  ' 
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We  also  introduced  the  abbreviation  of  a  vertical  bar  to  indicate  that  the  derivative  of  the 
overlap  matrix  is  to  be  taken  with  respect  to  the  R  dependence  of  one  side  only.  The 
derivatives  with  respect  to  nuclear  momentum  are  given  by  an  identical  expression  with 
the  gradients  V^.  Therefore  we  do  not  show  them  in  the  following. 

For  the  second  derivatives  of  the  overlap,  we  need  the  derivatives  of  the  inverse  of 
M.  From  the  relation 

XX~l  =  I  (A.10) 

for  any  matrix  X  depending  on  a  variable  x,  we  find 

^-  =  -X-^X-\  (AM) 
Ox  ox 

Therefore  the  derivatives  of  the  inverse  of  M  are 
and 


The  second  derivatives  of  the  overlap  are  then  given  by 

(A.14) 

-((r+.t.)-.t)  ((r+A)-.t) 
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and 


^ln(z|z) 
dz 


as  well  as 

<^ln(z|z)  _  ^ 
dz 


(A.15) 


The  derivatives  with  respect  to  nuclear  position  follow  in  the  same  way 


_d 

dz 


where  we  used  the  relation 

-i  / 


(A.16) 


-Vn\n(z',#\z,R)  l,=z,R,=R=(0    DVj  A|(/;)(/-  +  ,t2)-1 

-,(r  +  ,t2)-I(/.    2t)Vj?  A|(/2,)(/-  +  2t,)"1  (A.17) 


z^'  +  A)     =  (7°W)    z,  (A.18) 
which  follows  immediately  from 

(l0  +  zz^z  =  z(l*  +  z*zy  (A.19) 
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The  derivatives  with  respect  to  the  nuclear  coordinates  can  be  expressed  as 
=  Tr[(/-  +  ^)"V  ,t)V^A|(72') 


which  is  a  Hermitian  matrix. 


The  kernel  for  the  one-particle  density  matrix  or  1-matrix  is  defined  as 


Yix{z\z)  =  {z\z)-\z\b\b]\z) 


The  1-matrix  has  the  familiar  block  form 

r 

For  the  occupied  block  we  find  the  following  expression 

1-1 


/  P  T'\ 


r;t  ={z\z)~l(vae\ 


N    /  K  . 

/=1  \  m=N+l 


N    /  K  \ 

[bk  +    E  h°^mk)\vac) 

k=l  V  m=N+l  / 

=  (z\z)-1(-),+3det((Slk+    £  *ml*mk)¥i,ht)) 
\  m=N+l  ) 


=  det  (M)-1Minor(M)tJ 

=((^n. 


(A.21) 


(A.22) 


(A.23) 
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For  the  unoccupied  block  we  find 


rjt  =  {z\z)-l(vac\\(ztl-(b?+  £ 

A  V  m=N+l 


b°Jzml)b: 


N 


K 


II  K  +  E  h°m^l 
1=2  \  m=N+l 


AW  A-  \ 

II  K  +  E  ^mtlbac). 
k=2  \  m=A+l  / 


After  moving  the  annihilators  further  through,  we  obtain 


N  N 


/=1  *=1 


det  (  (6pq  +  2J  z*mpzmq)p^i^k 

V  m=l  > 


A 


=  det  (M)"1  ^  z^MinortM)^ 
Jt=i 


so  that 


(A.24) 


(A.25) 


(A.26) 


The  off-diagonal  block  is  obtained  with  the  help  of  the  first  derivatives  of  the  logarithm 
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of  overlap  kernel  from 

T'}l=(z\z)-\z\b?b)\z) 

=^Jln<*l*>  (A.27) 

Therefore  the  1-matrix  is  given  by 

r(z\z)=(rz^(r  +  z<zy\r  *t) 

(A.28) 

r(z'*,/2\P\z,i2,P)  =  (^(a'  +  AV'  +  AW^aVJ'V    *f ) 
Note  that  this  is  the  projector  onto  the  (non-orthonormal)  occupied  orbitals  [121]. 

The  two-particle  density  matrix  or  2-matrix  for  a  single  determinantal  wave  function 
is  a  simple  function  of  the  l-matrix[122] 

{zlz^iz^bihlz)  =  TkiTh  -  TkjTu.  (A.29) 


With  the  total  molecular  Hamiltonian  written  in  second  quantization  in  the  orthonor- 
mal  molecular  basis  as 

(A.30) 

K  K 
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we  find,  with  the  above  kernels  (A.23-A.27)  and  (A.29)  for  the  densities,  the  following 
expression  for  the  energy  of  the  coherent  state 

E(z*,z)  =£(°)  +  Tr^1))  +  ±Tt(vtW) 

=£(°)  +  Tv(hT)  +  iTr(Tr(Vat;o6r)ar). 


jTr(Tr(KaMar)ar)t 


(A.31) 


=E(°)  +  Tv(hT)  +  l-Tr(Tv(Vab.>abT)  T). 


2  /a  >b' 

Here  h  andV  are  the  one-electron  and  and  symmetrized  two-electron  integrals  in  the  spin- 
orbital  molecular  basis 

Viy,kl  =(ij\\kl) 

=(ij\kl)  -  (ij\lk)  (A.32) 
=(ik\jl)  -  (il\jk), 


with 


(til*/)  =  J 


^*(r1)0*(r2)V'Jt(n)V'/(''2)J3  ,3 

 ^  :  d  nd  r2.  (A.33) 

H  -r2 


When  straight,  i.e.  not  antisymmetrized,  two-electron  integrals  are  used  the  factor  of  one 
quarter  needs  to  be  replaced  by  one  half.  And 

e«»  =  rA  +  r  4*£ii_  (a.34> 

2Mfc    &      -  Ri\ 

is  the  classical  nuclear  kinetic  and  nuclear  repulsion  energy.  When  the  nuclei  are  treated 
quantum  mechanically  this  term  has  to  be  replaced  with  the  expectation  value  in  the 
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appropriate  nuclear  coherent  state  and  the  one-  and  two-electron  integrals  have  to  be 
replaced  similarly. 

The  one-electron  integrals  form  a  Hermitian  matrix.  The  two-electron  integrals  have 
the  following  symmetries: 

(ij\kl)  =(ji\lk) 

(A.35) 

=(ki\ijy  =  (ik\jiy. 

When  real  atomic  orbitals  are  used,  there  are  additional  symmetries  so  that  we  have  the 
following: 

(tj|*Q  =(ji\lk) 

=(kl\ij)  =  (lk\ji) 

(A.36) 

=(kj\il)  =  (jk\li) 
=(il\kj)  =  (li\jk). 

In  order  to  obtain  an  explicit  expression  for  the  energy  in  terms  of  the  coherent  state 
parameters  z  we  write 

=EW  +      kijTji  KWO  -  (iWWkiTij  (A.37) 

=EW  +  EW(z*,z)  +  EM{z*,z). 
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The  one-electron  energy  is  then 


&l\z\ z)  =Tv(hT)  =  Tr  (h( )  (/*  +       ~\  I'    z<  )j 
=Tr(AT)  +  2Re(Tr(/iT't))  +  Tr{h°T°) 
=TtU*(r  +  «t«)"1)  +  2Re^Tr^'z(/'  +  z^)"* 


(A.38) 


For  the  two-electron  energy  we  find 

*V,*)  =^Tr(Tr(v;6;a6r)ar)fc  =  ^W«r« 

=Wiv(v(lMi(^(/'  +  *t,)"1(i«  ,t)^  (A.39) 

(^)(r  +  .t,)>  .t,); 

Worked  out  in  blocks  this  becomes 

£(  V,  z)  J-  E  Viwnn  +  E  VwVN  +  \  E  KwHrirj 

+  E  +  Re(E  (A-4°) 

+  2Re(£  V^li,.rj.)  +  2Re(£  V^r^), 
where  we  have  used  the  symmetries  (A.35)  of  the  two-electron  integrals  to  reduce  the 
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16  terms  to  7.  We  conclude 


E(z*,z)  +  Tr(h(^rz^(r  +  z*z)  \p  z*)) 

+  W<&(v-g-^Vr  +  ***)"1(i*    z*))a  (A.41) 

(7;)(r+A)_1(7, 

Using  the  derivatives  (A.  12  and  A.  13)  of  the  inverse  of  M,  we  can  compute  the 
derivatives  with  respect  to  i  of  the  one-electron  energy  (A.38) 

^r=l4('(';)(-*-)-»-  ••>). 


K 

+  ' 


((o  i-^^JCr  +  .t,)"1) 


(A.42) 


+ 

Using  the  same  relation  as  in  Eq.  (A.  17),  we  find 

dE^(z\z) 


((,.,«.)-(-*   n^)(r  +  ^-')t.  (A.43, 


Applying  the  same  calculation  to  the  double  trace  and  using  the  symmetries  (A.35 
and  A.36)  of  the  antisymmetrized  two-electron  integrals,  we  obtain  the  derivatives  with 
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respect  to  z  of  the  two-electron  term  (A.39).  With  the  one-electron  term  this  gives 


(A  +  Tr(v^(^)(/,  +  *t*)"1(/-  )  (A.44) 

C)e+A)1- 

We  recognize  the  Fock  operator 

F  =  /t  +  Tr(VaM6r)a 

(A.45) 

=  ft  +  Tr([(oa|66)  -  (afc|6a)]r) 

in  this  expression. 

The  derivatives  with  respect  to  the  nuclear  coordinates  R  and  P  are  obtained  by 
taking  the  derivatives  of  the  one-  and  two-electron  integrals  h  and  V  and  of  the  overlap 
matrix  A.  The  derivatives  with  respect  to  z  are  the  same  whether  we  consider  orbitals  that 
depend  on  nuclear  positions  or  not,  and  for  that  reason  we  have  not  shown  the  dependence 
of  E  on  the  nuclear  coordinates  and  momenta.  The  expressions  for  the  derivatives  are 
Eq.  (A.37-A.41)  with  the  integrals  replaced  by  the  corresponding  derivative-integrals, 
plus  the  reorthonormalization  terms  involving  the  derivatives  of  the  overlap  matrix.  All 
derivatives  with  respect  to  nuclear  coordinates  are  understood  to  be  derivatives  with 


(A.46) 
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respect  to  R',  F  which  ,  afterwards  are  set  equal  to  R  and  P,  respectively,  i.e. 

!  ZkZte2(Rk-Mi) 
V*  E(z*,R,P,z,R,P)  =  -->   -  v   -  ,  ' 

+  Tr(v^ftr)-Tr(vAArAr) 
+  5Tr(Tr(v*.v'^r)/)J 

-Tr^VjArV,,,,^^. 

When  the  nuclei  are  treated  classically,  both  the  integrals  and  their  derivatives  with  respect 
to  R  are  given  by  available  programs.  The  derivatives  with  respect  to  P  are  similarly 

— ♦ 

VpE(z\R,P,z,R,P)  =^-  +  Tv{yphT)-Tv(ypAThT) 

+  iTr^Tr(vAVa6;afcr)ar^  (A.47) 

-Tr(Tr(vAArv;6;a6r)ar)6. 

The  terms  after  the  first  are  only  present  when  travelling  molecular  orbitals  are  used. 
Currently  no  programs  are  available  to  give  such  integrals  or  their  derivatives.  For  more 
advanced  nuclear  descriptions,  the  integral  evaluation  requires  a  more  detailed  treatment 
and  the  derivatives  get  correspondingly  more  involved. 


APPENDIX  B 
DETAILS  OF  DERIVATION  IN  ATOMIC  BASIS 

From  Eq.  (B.l)  for  the  coefficients  v  of  the  dynamic  orbitals  orthogonal  to  the 
occupied  ones,  we  define  the  following  matrices 

A»(*)  =  (J«   *t)^*   ^VQ=A*  +  ^At  +  A'«  +  «tA°z  (B.l) 

A  being  the  atomic  basis  overlap  matrix  and 

\°(v)  =  (v    P)(£    ^j^j=A°  +  ,A'  +  AV  +  VAV  (B.2) 

which  naturally  occur  in  many  expressions. 

From  the  expression  of  the  parameters  v  as  a  function  of  z  the  following  derivatives 
are  calculated 

Czmn  V  /  nj 

+  ((a'*  +  A°z)  (A"  +  A'*)"1  A')  +  A'z)-1)  (B.3) 

=  -(A°  +  ^)im((A'  +  A',)-1)^ 

^=((A.  +  ^)-)j>(At(A.  +  ^)-I(A'  +  ^))my 

-((a'  +  ^A*)"1)   A^  (B.4) 

=  -((a-+^)-')  (a'+aV)^. 
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As  in  the  previous  appendix,  the  derivatives  of  the  logarithm  of  the  overlap 


5  =  det{(7-  *f)A('*)} 
=  det(A*  +  A'z  +  2fA'f  +  z^A°z) 
det(A*(z)) 

have  to  be  computed.  Repeating  the  calculations  in  (A.2-A.8),  we  obtain 

~-ln(zV>U'=«  =  E  (A'  +  zlA°)tk  (A*  +  A'2  +  ztA'f  +  z^°z)  -i 
-(A-'(A'  +  ,tA-)) 


(B.5) 


(B.6) 


^In^'l^l^,  =  T  ( A'*  +  A°z)    (A*  +  A'z  +  zW  +  ztA0z) 
=  ((A'.  +  A°,)A-)ti 


and 


Valn^i^P'I^P) 


z'=z,R'=R,P'=P 


=  TV[(A'  +  A'z  +  2fA"  +  2fA0z) 


-1 


(B.7) 


(B.8) 


K A'l  +  *A  A'l*  +  **Vj  A"l  +  *fVA A°l*)] M 
and  an  identical  expression  for  the  derivative  with  respect  to  nuclear  momentum. 
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The  derivatives  of  the  inverse  of  overlap  matrix  of  the  occupied  dynamic  orbitals 
are  similar  to  (A.  12  and  A.  13) 


d 


dzkj 


and 


A*  +  A'z  +    A'f  +  z^A°z\  J 

/  mn 

--(A^r^A'  +  jftA-))  (B.9) 


^-((A'  +  A'^  +  ^A't  +  ztA^) 

=  -fA*(^)-1)  (B.10) 

V  /  mi 

((At  +  A-»)A-W-') 


The  second  derivatives  are  then 


(A'(.)->(A'  +  .tA«)). 


and 


=  (A0  -  (A't  +  A^A'^)"1^'  +  jJA-JMA-M-1) 


(B.ll) 


=  £  j^A'W^A'  +  .tA-^i  +  (A^)"1)^ 

,=1  ^/m  (B  12) 
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Using  the  virtual  dynamic  orbitals  v  given  by  (3.46)  and  the  relation  /  -  Pocc  =  Pvirt 
expressed  in  terms  of  the  atomic  orbitals,  then  J  —  Pocc  is 

(A9  A'\ 

~\A*   A0  J' 

(  (A*  +  A'z)A*{z)-1  (A'  +  ^A't)     (A*  +  A'*)A#(*)_1  (A'  +  z^A°)  \ 

'  V  (A't  +  A°z)A'(z)-1  (A*  +  zf A't)    (A*  +  A°z)A*(z)~1  (A'  +    A0)  J 

(B.13) 

and  the  virtual  projector  is 

(fvt£o)A»-VA-  +  A't  + 

/  (AV  +  A')A°(u)_1  (vA*  +  A'f)     (AV  +  A')A0(u)_1(uA'  +  A0) 
=  V  (A^v*  +  A0)A0(v)_1  («A*  +  A't)    (A'V  +  A^A^v)"1^'  +  A0) 

(B.14) 

Then,  from  the  equality  of  these  two  expressions  we  can  simplify  Eq.  (B.12)  to 


and 


g=((.V  +  i')A»-'K  +  A.))n 


g8--((^  +  A-.)A^)ta 


(B.15) 


((A^+A'ijA'W"1)  . 


(B.16) 


100 


The  mixed  derivatives  with  the  nuclear  coordinates  are 


-  (A"  +  tfz)\\z)-\r  z^V^Al^y^z)-1 

The  double  derivative  with  respect  to  the  nuclear  coordinates  is  given  by 

V^V^ln^i^,/?)^^ 


=  Tr 


A*(z)-\r  **)VAV4A|I 


(B.17) 


(B.18) 


The  1-matrix  (Eqs.  (A.23)-(A.28))  becomes  in  the  atomic  basis 


rat=A^/2*^(A,  +  A'z  +  2tA't  +  ztA°z)  \l*  zf)A 


=ArA. 


(B.19) 


Notice  the  overlap  matrices  on  both  sides  as  is  typical  for  projection  operators  in  a 
nonorthogonal  basis.  T  is  the  operator  form  of  the  1-matrix,  whereas  Tat  is  the  matrix 
element  form  It  is  typical  in  nonorthogonal  bases  that  every  linear  operator  has  two 
matrices  representing  it.  These  are  the  covariant  and  contravariant  forms. 
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Then  the  energy  (A.41)  of  the  coherent  state  becomes 


(B.20) 


E{z\  z)  =f?(°)  +  Tv(hT)  +  ^Tr(Tr(K6;air)ar)6 

+  iTr(Tr(va6;a6(/;)A-(z)-1(/-  ^ 
(7;)a-(2)-(/.  .t^ 

which  is  only  slightly  more  complex  than  in  the  orthonormal  basis! 


For  the  derivatives  of  the  energy,  only  the  derivative  of  the  one-electron  part  is 
shown  in  detail,  the  two-electron  part  is  similar.  Using  (B.9  and  B.  10),  the  derivatives 
of  (B.20)  give 


dzii 


dEiy  >*)  (At  +  A«z)  (A-  +  A'z  +  ztA't  +  .tAO,)-1 

(A'  +  zU't    ti  +  2ih°)  +  {h't    h°)^j  (B.21) 


(^(a'  +  A's  +  z^  +  ^z)  ^  . 


Using  the  Eqs.  (B.13)  and  (B.14)  with  the  rightmost  A  replaced  by  the  one-electron 
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integrals  h 


h°  -  (A't  +  A°z)  A'(z)-1  (h'  +  ( Pz  )  A'(z)-1^ 

=  (((A'V  +  A0)  A'ivy1  (vh'  +  A't) 

(AV  +  A0)  A'iv)-1  (vti  +  >>0))  ( ^  )  A*(^)-1) 
=  ^(A'V  +  A0)A°(l;)-1(t;    P^'^A'iz)-1^  . 
With  a  similar  derivation  for  the  two-electron  part  we  find 


(B.22) 


^f)=((AV+A-)A»-'(»  /•) 

c>>-), 

with  v  given  by  Eq.  (3.46).  We  again  recognize  the  Fock  operator 

F  =  h  +  Tr(Vab.tabT)a 

=  h  +  Tr([(aa\bb)  -  (ab\ba)]T)a. 


(B.24) 
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The  derivatives  with  respect  to  the  nuclear  coordinates  are 

V.  E(z*,B!,P',z,R,P)  =-^Y  ^-=t  - 

Rk    V  1  R'=R,P'=P        2^  \R,-R,\3 


and 


\Rk  -  Ri\2 


+  Tr(vAhr)-Tt(vjUArhr) 

+  H*(vAW).r)| 
-Tr^Tr(v^Arva6;a6r)ar^ 


VAE  =  K  +  Tr(vAftr)  -  »(vaat*t) 


(B.25) 


+  ilr^Tr(vAVfl6;air)ar^  (B.26) 
-Ir^Tr^Ar^r)^  . 
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Finally  he  did  an  about  face  and  decided  to  work  on  the  very  small,  but  not  the  tiny. 

It  is  of  small  consequence,  but  interesting  to  note,  that  at  the  age  of  7  he  asked  his 
mother  what  jobs  there  were  which  did  NOT  use  mathematics,  since  he  simply  hated  the 
subject.  That  was  to  change  only  after  he  was  introduced  to  algebra. 
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